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ABSTRACT 


Current research at Rensselaer Is generating fun- 
damental engineering design techniques and concepts for the 
chromatographic separator of a chemical analysis system for 
an unmanned, martian roving vechlcle. Previously developed 
mathematical models of the gas chromatograph were inadequate 
' for predicting peak height and spreading for som^ experimen- 
tal conditions and chemical systems* Consequently, a new 
chromatograph model Is developed which Incorporates previously 
neglected transport mechanisms* A closed— form analytical 
solution to the model Is not available so the numerical tech- 
nique of Orthogonal Collocation is studied. To establish the 
utility of the method, three models of Increasing complexity 
are considered, the latter two being limiting cases of the 
derived model; 1) a simple, diffusion-convection model; 2) 
a rate of adsorption limited, inter-lntrapartlcle model; and 
3) an Inter-lntrapartlcle model with negligible mass transfer 
resistance. The first model involves one dependent variable 
and one spatial dimension; the second, two dependent variables 
and one spatial dimension; and the third, three dependent var- 
iables and two spatial dimensions. The orthogonal collocation 
treatment reduces the models to sets of ordinary differential 
equations which are integrated using the Bulirsch-Stoer ex- 
trapolation technique. 

Simulations with the first model using actual chro- 
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matographlc Input pulse data show the collocation procedure 
to accurately represent system behavior* Large Peclet numbers 
usually observed In practical chromatographic columns require 
a higher degree of approximation than low values* In general, 
15 collocation points suffice* Similar results are obtained 
from a study of the second model which Involves two coupled 
partial differential equations* The model Is successfully 
solved numerically, although computation time becomes exces- 
sive* The investigation Is concluded with a preliminary 
study of the third model which involves three coupled partial 
differential equations. Estimated computational times based 
upon partial simulations of this model show complete numerical 
solution within available computer capabilities and financial 
constraints to be unfeasible. It. Is therefore concluded that 
if orthogonal collocation is to be applied successfully to 
pulsed, distributed systems of the chromatograph within com- 
puter constraints, further research on the different charac- 
teristics of the orthogonal functions and the formulation 
the trial function must be undertaken. 
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PART 1 




INTRODUCTION AND SUMMARY 

The synthesis and analysis of mathematical models 
of the gas chromatograph is one subtask of a group effort de- 
signed to define fundamental design criteria for an optimal 
design of a combination gas chromatograph - mass spectrometer 

V 

system which is to be part of an unmanned mission to Mars* 

The task which must be performed by this part of a Martian 
Roving Vechicle Is the analysis of samples to determine the 
existence of organic matter and living organisms on the 
Martian surface. The analysis will involve the subjection of 
gaseous,, liquid and solid samples to biological and chemical 
reactions, with subsequent product separation and identifica- 
tion using the gas chromatograph - mass spectrometer system. 

The chromatograph may be looked upon as a separating 
device where, the phenomenon of adsorption-desorption Is util- 
ized. Owing to the different characteristics of various chem- 
icals, each species will adsorb and desorb at different rates 
when exposed to a packed bed of granuaar particles with or with, 
out a liquid substrate. Because of the unique behavior of 
each chemical, a multicomponent sample may be injected Into a 
chromatograph and elute as separate waves of specific chemical 
species. 

prior to this Investigation, chromatograph models 
have been formulated based on Interparticle transport mechan- 
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isms with simple adsorbed phase behavior assumed* These pre- 

I 

vious model formulations, which have all had closed-form, 
analytic,' time-domain solutions, have proven incapable of ade- 
quately predicting component behavior in all cases, Conse- 

f . 

quently, a new model has been derived which Includes both. 
Interparticle and intrap article transport mechanisms. The 
complexity of this new model precludes direct analytical sol- 
ution, and hence application of an appropriate numerical tech- 
nique is necessary to effect time-domain solution. Prior to 
any time -domain analysis, the model is analysed In the Laplace 
transform domain using the method of moments. The first three 
moments of the Impulse response of the model are derived. 

Using actual input data, predictions for the first three moments 
of the output data are made and are compared with actual 
chromatographic data and predictions of a simpler, interp article 
model. The results indicate that the new model is more ca- 
pable of predicting the momenta, of the actual data. 

Because the mathematical complexity of the new model 
prohibits a direct, closed-form analytic expression for a res- 
ponse, appropriate numerical techniques applicable to the equa- 
tions of the new model (and future models which may involve 
nonlinear terms) must be used to allow direct comparlsions 
between prediction and experiment. For the systems of equations 
encountered in chromatograph modeling, numerical techniques 
require a finite terminal boundary condition as opposed to an 
infinite column boundary condition usually used in deriving 
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analytical solutions to simpler chromatograph models. As a 
prelude to numerical technique considerations, a simple, tran- 
sient, diffusion-convection, mass transfer- equation is analysed 
and criteria are developed wherein a finite terminal boundary 
condition can be applied to yield Infinite column behavior at 
the bed outlet. 

Subsequent to the development of a complex chromat- 
ograph model and the realization that model simulation and ver- 
ification require a numerical technique, attention is directed 
to the study of -Orthogonal Collocation as a technique suit- 
able for routine analysis of complex chromatograph models* The 
motivation for conducting this investigation is several-fold: 
computational limitations of the widely used Finite Difference 
method, successful use of orthogonal collocation to solve cer- 
tain chemical reaction engineering problems, and the lack of 
documentation for the application of orthogonal collocation to 
pulsed, distributed systems such as the chromatograph system. 

The general theory and computationpl aspects of orthogonal 
collocation are reviewed and discussed. One of the steps in 
the application of orthogonal collocation Involves the inversion 
of a matrix. Previously documented developments have given 
formulations where the matrix to be Inverted becomes increas- 
ing ill-conditioned with Increasing size and may, due to com- 
puter precision limitations, prove non- Inver table. Hereto- 
fore, this has not been documented. An alternative develop- 
ment, theoretically equivalent, is presented which is shown 









to successfully eliminate this problem to a high degree. 

In order to establish whether orthogonal collocation 
is a technique worthy of exploitation In the analysis of chro- 
matographic systems, three models of Increasing complexity 
are solved using the methodj 

Xm A simple, transient, diffusion-convection mass 
transfer problem. 

2, A rato-of -adsorption-limited, Inter-intraparticle 
problem - a limiting case of the derived model. 

3» An inter-intraparticle adsorption problem with 
negligible mass transfer resistance between the 
interpartlcle and intraparticle regions - another 
limiting case of the derived model. 

Por each model, the orthogonal collocation treatment reduces 
the partial differential equatlon(s) to a set of ordinary 
differential equations. 

The first model is studied because It possesses re- 
sponse characteristics found In more complex models, possess s 
an analytic solution for direct comparlsions with numerical 
results, and establishes guidelines for more complex models 
to be considered. Prior to solution of the resultant set of 
ordinary differential equations, an eigenanalysis is made of 
the differential equation set. This set or the resultant dis- 
cretization of the distributed system Is stable for axial Peclet 
numbers from 1 to 10000 and approximation orders of 3 to 21. 

This model Is solved for cases of rectangular and actual system 
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Input data pulses. The effects of high (10000) and low (32) 

i 

values of axial Peclet number are studied to determine the de- 

! 

j 

gree of approximation required for good representation of the 
exact system response. The high Peclet numbe^ (10000), char- 
acteristic of the magnitudes encountered in actual chromato- 
graphic system data, requires a higher degree of approximation 
than the low value. The sharpness of the forcing function also 
affects the numerical results; t a higher order of approx- 

imation is required for very sharp Input pulses. For the siroother 
actual input data and the high Peclet number, a fifteenth or- 
der approximation is adequate. The set of ordinary differential 
equations is Integrated using the Bullrsch and Stoer extrap- 
olatlon technique. For this type of problem, this method is 
more efficient (for equivalent error tolerance) than the more 
well-knovjn Euler, fourth order Runge-Kutta and Hamming Predictor- 
Corrector techniques. Consequently, the extrapolation treat- 
ment is used exclusively for Integration of the sets of ordin- 
ary differ. ntlal equations that result from the application 
of orthogonal collocation to the problems considered in this 
Investigation. 

Following a study of the simple model, orthogonal 
collocation is applied to solve the second model given above. 

This problem possesses an analytical solution which is used 
for comparision with the different degrees of approximation 
considered. The system parameters which appear in this model 
correspond to parameters encountered in actual chromatographic 



6 


system experimental work# This is Important because the 
parameter choice; e*£. • Peclet number. Is dictated by actual 
experiment rather than convenience* This problem is more 
complex than the simple problem in that two coupled partial 
differential equations are treated using the orthogonal collo- 
cation method* As a consequence, a higher degree of approxi- 
mation is necessary and the constraints of excessive computer 
time and suitable computer hardware availability come to the 

forefront of the investigation, 

The investigation is terminated with the application 
of orthogonal collocation to the third model listed above# 

This model has no direct, analytic solution available. Hence, 
the strategy is to successively apply orthogonal collocation 
with increasing orders of approximation until a convergent 
response Is realized* Unlike the previous two problems where 
only one spatial domain is discretized, the interparticle, 
this problem requires orthogonal collocation discretizations 
for both the interparticle and intraparticle regions. The 
problem involves the solution of three coupled partial differ- 
ential equations. Again, actual chromatographic system param- 
eter: values are used# To effect the above strategy of succes- 
sive simulations with increasing order of approximation, sever- 
al cases are studied for short computer run times# These ' 
times are extrapolated to give estimates of computer require- 
ments necessary to complete the analysis# These extrapolations 
indicate that within available computer hardware capabilities 



and financial constraints, thorough analysis of this problem 
ts not feasible using orthogonal collocation. This does not 
rule out the utility of the theory of orthogonal collocation 
as a technique but points out a problem where innovation and 
further study may be necessary for the realization of a prac- 


tical solution 



PART 2 


CHROMATOGRAPH SYSTEM MODELING 

A# Chroma top;r a ph Modeling Background 

One area of the overall gas chromatograph systems 
study has been the mathematical modeling of the chromatograph 
system* Work In the area has been carried out by several 
investigators (Silva, 1968} Voytus, 19^9; Taylor, 1970; Keba 
and Woodrow, 1972)* A course has been pursued wherein succes- 
sively more complex models have been considered. These 
models have all yielded analytical expressions from which a 
simulated chromatogram could be computed directly, Compar- 
ision of predicted system behavior with actual system data 
has directed modeling efforts to consider more adequate and 
hence more complicated models, 

“Prior to this investigation, the most complex model 
proposed for the chromatograph system was based on an Inter- 
partlcle phase mass balance and an adsorbed phase mass balance. 
Several transport mechanisms were included; axial diffusion, 
convection, and mass transfer between the interpatlcle and 
adsorbed phases. A linear isotherm was used to describe the 
adsorption kinetics. This model has been studied and compared 
(Keba and Woodrow, 1972) for the cases of finite rates of 
mass transfer to the adsorbed phase (nonequilibrium adsorption) 
and infinitely high rates of mass transfer to the adsorbed 
phase (equilibrium adsorption). In both cases, simulations 
using the models failed to predict the degree of dispersion 
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exhibited by many of the experimental data. It was concluded 
that additional transport mechanisms, e.g. , ihtrapartlcle 
diffusion, may be contributing appreciably' to | the overall ad- 
sorption-desorption process. Hence, further model develop- 
ment and analysis was indicated, 

B. Development of the Inter-Intraparticle Adsorption 
Model 

Previously, the intraparticle region of the chromat- 
ograph packing material has been modeled as being nonexistent 
or as a region where the transport processes occur at such a 
rapid so as not to significantly affect the dynamic behavior 
of the system. It is the purpose of this section to refor- 
mulate the chromatograph system model by including the trans- 
port process which are presumed most likely to affect the 
dynamics of the adsorption-desorption process within the 
chromatograph, packing material. 

Figure 1 presents graphically the transport pro- 
cesses to ''e modeled. The sample to be separated is injected 
into a relatively inert carrier gas, e.g., helium. As this 
slug of sample is transported down the chromatograph by the 
carrier gas, the various species diffuse, adsorb, and desorb. 
Diffusion of the chemicals in the direction of the carrier 
gas flow in the interparticle region is represented by the 
dimensionless parameter, which is determined by the 

system fluid mechanics. Mass transport from the Interparticle 
region to the intraparticle region is represented by a 



FIGURE 1 
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dimensionless parameter, N^qq, which is essentially deter- 
mined by the system fluid mechanics. Diffusipn in the Intra- 

i 

particle region is represented by a dimensionless parameter, 

> ! 

Pey^, which is in part determined by the properties of the 
particle packing. The rate of adsorption within the particle 
is characterized by the dimensionless parameter, Ngy* Adsorp- 
tion-desorption within the particle is represented by mRi, a 
thermodynamic parameter peculiar to each species. This param- 
eter contains an equilibrium constant, m, and the quantity 
Rj. Rj is the ratio of moles of fluid within the particle to 
the moles of adsorptive sites within the particle. The quan- 
tity Ri Is directly related to the quantity Rq where Rq Is 
the ratio of moles of fluid within the total bed to the moles 
of adsorptive sites within the bed. The relationship between 
these quantities is: 

Ri = (€/(i-6)ya) Ho (i)» 

The reason for noting this relationship is that the parameter, 
mR() has been noted in previous models and the above relation- 
ship serves as a unifying concept for the new model formul- 
ation which follows. 

With the above concepts in mind, the following set 
of dimensionless equations has been derlved^^* based on the 
assumptions which follow: 

See Part Nomenclature, for definition of terms, 

** See Appendix A for derivation. 
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An Intraparticle phase mass balance ;! 


An adsorbed nhase mass balance: 



Ng„(yi - y*) 


( 2 ) 


( 3 ) 


( 4 ) 


A thermodynamic relationship between the intraparticle 
and adsorbed phases : 


yj = m Xa 


( 5 ) 


The above equations are valid under the following assumptions: 

1, The column is isothermal. 

2. The carrier gas velocity profile is flat. 

3* The axial diffusion coefficient is a composite 
factor which may or may not have a turbulent 
component. 

The gas composition Is: approximately constant 
in the radial direction at a given axial posi- 
tion. The concentration gradient occurs in a 
thin boundary layer at the inter-intrapartlcle 
Interface. 

5 * The gas composition within the particle is ap- 
proximately constant in the angular direction 





at a given radial position; the concentration 
gradient occurs only in a thin boundary layer 
near the adsorbent surface. 

6. The adsorbent layer is so thin that there is 
no diffuslonal resistance within the layer in 
the direction normal to the surface. 

7 . The diffusivity in the adsorbent layer is so 
small that there is no diffusion in the direction 
parallel to the surface in the intraparticle 
radial direction. 

8. The net rate of adsorption for the carrier gas 
is negligible. 

9* Only one component is adsorbed and its gas phase 
composition as a mole fraction is small compared 
to unity. 

10. ' The carrier gas behaves as an ideal gas. 

An applicable set of boundary and Initial conditions 


are as follows : 

Initial Conditions ; 

y (z, 0) = 0 (6) 

yi 6 } = 0 (7) 

Xa 0) = 0 (8) 

Boundary Conditions ; 

•y(0, e) = (9) 

[(^/6 )apL(L/R)/Pe^] when /z=l (10) 



hy^/6rL = 0 ; /v= 0 (11) 

11m y(z, G) = finite (12) 

These conditions reflect a sample-free column at zero time, 
a sample injected as an impulse, mass transfer between the 
Interparticle and intraparticle regions, no concentration 
gradient at the center of the column packing, and no end 
effects at the column exit. 

For the systems under consideration it has been 
shoim by Keba and Woodrow (1972) that inclusion of the par- 
ameter NtoG Is of minor importance. If one were to consider 
the case of infinite rates of mass transfer, l.e., 
the coupling condition given by equation (10) would be re- 
placed by 

yi (z» 1. e) = y (z, e) (13) 

Thus, a model in the form of a set of coupled, 
partial differential equations is proposed. Prior to con- 
sideration of the time domain solution of the equations, a 
moment analysis can be made to ascertain the predictive ca- 
pabilities of the proposed model. This analysis is the subject 
of the next part of this investigation. 
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MOMENT ANALYSIS OF THE INTER- INTRAPARTICLE 
ADSORPTION MODEL 

A, Theory and Backfcround 
An analysis of a proposed model can be made prior 
to determination of the model’s time-domain solution to 
yield the gross characteristics of the impulse response of 
the model. In addition, because of the poor predictions of 
previous models (Keba and Woodrov/, 1972) with respect to 
chromatogram spreading, it is desirable to know the nature of 
the response of the proposed model for the pulse-type functions 
which are the sample injections seen in experimental work. 

The nature of the response can be characterized by statis- 
tical quantities known as moments which may be obtained v/ith- 
out knovfledge of the time -domain model solution. The moments 
may be derived directly from the Laplace domain solution of 
the model. The following development will indicate how the 
moments of a model are obtained and hov; the analysis can be 
extended to give the moments of systems forced by general 
pulse-type Inputs. ' 

The Impulse response of the chromatogram may be 
viewed as the residence time frequency distribution (Douglas, 
1972). This quantity resembles the probability distribution 
function which appears in statistical analysis. The moments 
of the distribution function about the time origin are defined 
by the following: 
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-«n = /e" 

o 

where 

f(0) = the distribution function being analysed. 

The denominator of equation (14) is the area under the function. 
The relationship of the moments about the origin to the 
Laplace transform is developed in Appendix B. The result is; 




(15) 


( 16 ) 


Interest also centers on the moments about the 
first absolute moment or mean, Mathematically these are 

defined by: 

M 

These moments about the mean A/, are directly related to the 
moments about the origin. The relationships are obtained by 
formal expansion of equation (1?). Appendix B gives the rel- 
ationship for n=2 and n=3* For n=2, the moment about the 
mean is exactly the variance of the response. For n=3, the 
moment about the mean is related to the skew of the response. 
One can use the preceding to develop equations rela- 


^ 

n = / f(e)de / J f(e)de ; n>; 


( 1 ?) 


ting the moments of system responses for arbitrary pulse-type 
forcing functions (see Appendis B for details). That is, 
given the system input data (the moments of which can be. 



1 ? 


computed from equations (14) and (1?)) and the system trans- 
fer function (the Laplace transform of the Impulse response), 
the moments of the system response may be determined and 


compared with the moments of the actual output data. Referring 


to the block diagram in Figure 2, the results are; 

Ay = Ax ' Aq 

^ly = ^IX *** 

/^2Y = ^2X ^2G 

MjX “ ^3X ^3G 


(18) 

(19) 

( 20 ) 
( 21 ) 


Equation (18) states that the area under the output curve is 
the product of the area under the input curve and the impulse 
response curve. Equation (19) states that the mean of the 
output occurs at the sum of the mean of the input function 
and impulse response. Equation (20) states that the variance 
of the output is the sum of the variance of the input function 
and the variance of the impulse response. Equation (21) states 
that the third moment about the mean of the output is the sum 
of the third moments about the means of the input function 
and impulse response^ respectively. 

This technique can also be used for estimating sys- 
tem parameters. Douglas (19?2) uses an equation similar to 
equation (20) to estimate an axial Peclet number for a packed 
bed, Schneider and Smith (19<^8) apply moment analysis to es- 
timate adsorption equilibrium constants, rate constants, and 
intraparticle dlffusivitles for a chromatographic system mod- 
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X(S)= [x(0)] •, X(9) IS THE FORCING FUNCTION 

Y(S)= ^[yO)] i Y(0) IS THE SYSTEM RESPONSE 
G{S)= SYSTEM TRANSFER FUNCTION 


FIGURE 2 TYPICAL SYSTEM BLOCK DIAGRAM 
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eled similarly to that of Part 2, However, accurate param- 

i 

eter estimation using this method is limited by the accuracy 
of the data used for analysis. | 

I 

B. A-pplication of Moment Analysis to the 
Inter-Intrapartlcle Adsorption Model 
^The previous section outlined a method which can be 
\ used to analyse pulsed systems to determine the gross charac- 
teristics of the system response. This section will document 
an application of the concepts of moment analysis to the pro- 
posed model of Part 2. 

Consider the set of partial differential equations, 
boundary conditions, and initial conditions, equations ( 2 ) 
through (12). A Laplace transform domain solution for the 
Impulse response or transfer function was derived and appears 
in Figure 3; details appear -in Appendix C. 

Applying the definition given by equation (15) and 
using equation (1?), the moments A,, /(j, and are derived 
for the Impulse of the Inter-Intraparticle Adsorption Model, 
The results are presented in Figure 4; details of the manipul- 
ations appear in Appendix D. 

The parameters Peg, J^®A estimated 

a priori . The parameters mRQ and are not predictable a 
priori . Previous modeling analysis has estimated mHQ by a 
curve fitting process (Benoit, I 97 I). The estimation of 
will most likely Involve curve fitting also. 

An analysis was made using existing single component 
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y(i,s) 


where : 

K(s) 

X(s) 


ai(s) 




exp 



"X(s) ) + s 

b slnh(yax ) 

[(b-1 )sinh(>/ax ) + ^/ai cosh(v^)3 


-Nru mRx 
(s + NjiU niHx) 


%U ® 


,r*2 


(£) 



2 1 
Pe^ 


Particle porosity 
Bed void fraction 


Figure 3 Transfer Function for the -Inter- Intra particle 
Adsorption Model 



M.X I+I/idRq + (1 — 

M2 = 2(yt<i)2/PeE + 2-[(l-€)^/e] 

♦ m + l/niHi)2 [(R/L)2 PeA /15 + (l-O,^/£Nt0G] 

+ I/NeuCiJiRi ■ 

= 6M1M2 /Peg ■•■ ^ ^ [(l+l/inRj)/N2u(B®l)^3 

[(R/L)^ PeA /15 + (l-6)P/eNtoG] + (l+l/mRi)3 

L((l-e)/9/eNtoG)^+ 2(l.£)/e(R/L)2peA/15eMtoG 
- 23 (R/L)^Pg|/315] + 1/Kru (oEi)A 


Fi^re 4 Moments of the Impulse Response of the 
Inter-Intraparticle Adsorption Model 
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data. The parameters Peg, Pe^* and K^qG estimated using 

existing correlations. The values of ihRq which were estimated 
by Keba and V/oodrow (1972) using simpler models were used and 
the Ngu was varied. Tables 1 and 2 give results of this anal- 
ysis for acetone at 100 C and ethylene at 50 0, Both exper- 
iments. used Chromasorb 102 column packing, a porous material. 

In each case, the moments for the impulse response of the model 
were computed using the equations given in Figure 4. Use of 
system input data and equations ( 19 ) through (21) give predic- 
tions as a function of for the output moments. These pre- 
dicted values are compared with actual moments of the output 
data and with the predictions of the simpler, interparticle 
equilibrium adsorption model. Expressions for the moments of 
the simpler model were initially developed by Voytus (19^9 )• 

The results indicate that the proposed model can' 
more closely predict the characteristics of the output data 
than the simpler, interparticle model. The results indicate 
that a value of order of several hundred will give 

a predicted second moment very close to the second moment of 
the output data. This magnitude of is consistent with the 

values of which can be deduced from the independent research 
of Schneider and Smith (I 968 ). Tables 1 and 2 further indicate 
•that matching of the third moments would give different values 
of However, the use of third moments is not as reliable 

because data inaccuracies are further magnified in the analysis. 

It should be noted that if one accepts the value of 



TABLE 1 


MOMENT ANALYSIS AND PARAMETRIC 
STUDY - ACETONE lOO'C. 


mRo 

0.029 


, observed 
173.29 


M. 


( 1 ) 

1 , predicted 

I5O9 


M 


( 2 ) 

1 , predicted 



^RU ^2, observed predicted 


100 815.67 977.55 

200 , 

723.41 

400 

686.34 

800 

517. 80 

1600 

483.53 

3200 

466.40 

6400 

457.83 

12800 

453.55 

2S600 ' 

451.41 


^2, predicted 
'Ij-37.26 



^BU >^3, observed -^3, predicted ^3t pi'edicted 

166 25404.0 23192.0 19499.6 

200 

20454.8 


400- 

19745.3 


800 

19555.4 


1600 

. 19501.7 


3200 

19485.1 


6400 

19480.4 


12800 

19477.2 


25600 ' 

19476.2 1 



Peg = 8689. 

KtoG = *^8960. j 

(L/R)Vpca =328.2 

(1) Inter- Intraparticle Adsorption Model 

(2) Interparticie Equilibrium Adsorption Model 
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TABLE 2 


MOMENT ANALYSIS AND PARAMETRIC 
STUDY - ETHYLENE 50 C. 


mHo 

, observed 

.predicted 

(2) 

‘”1 , predicted 

0.194 

26.475 

257956 

23.719 


^RU -^2, observed ^2, predicted '^2, predicted 

100 7.024 13 . 2 S 3 o75^B 

200 

6.973 


400 

3.61? 


800 

2.240 


1600 

1.451 


3200 

1.056 


6400 

0.859 


12800 

0.760 


25600 

0.711 1 



^RU ^3 .observed '^3* Predicted /^3. predicted 

100 . 19.623 ~ 13.049 0.191 

200 

3.519 


400 

1.058 


800 

0.403 


1600 

0.219 


3200 

0.163 


6400 

0.144 


12800 

0.137 

■ 

25600 ! 

0.134 



Peg = 9744. 

NtOG = 79750. 

(L/R)^/Pej^ = 436.2 

(1) Inter-Intraparticle Adsorption Model 

(2) Interpartlcie Equilibrium Adsorption Model 
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Nj^U being on the order of several hundred for each case, 
all other parameters, excluding iiiRq, are of the same magni- 
tude* The key to the difference in the two component behaviors 
Is the parameter hiRq* 



PART 4 

TERf-IINAL BOUNDARY CONDITION ANALYSIS 


Mathematical modeling of chromatographic systems 
commonly require solutions to equations of the form: 

(l/Pe)(^V/^2^ ) « - Ra = dy/<36 (22) 

Application of analytical techniques to the above equation, 
when possible, commonly utilize the terminal boundary condition: 

11m y(z,e) = finite; G>0 (23) 

2-^00 

Use of the above boundary condition in analytical work yields 
a great deal of mathematical simplification. In addition, the 
use of this boundary condition is consistent with the theory 
which has been developed for prediction of the dispersion In 
packed beds; see, for example, Gunn (1969)* 

However, when numerical techniques must be applied 
to solve equation (22) or any other model which defies analyc- 
leal solution, the terminal boundary condition given by. equa- 
tion (23) must be replaced by a terminal boundary condition 
which is both computationally expedient and physically mean- 
ingful. A finite boundary condition v/hloh has found general 
usage in chemical reaction engineering problems (Danckwerts, 
1953 ). and (Hehner and Wilhelm# 195^) is: 

^y(i.e)/&z = 0; e>o ( 24 ) 
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Bastlan and Lapldus (195^) considered the case where R In 

equation (22) was an adsorption term. A llnelr relationship 

i 

I 

was assumed to describe the adsorption kinetics. For a step- 

I 

input and the conditions chosen, Bastlan and Lapldus showed 

j 

that finite column calculations, using equation (24) as a ter- 
minal boundary condition, closely approximated infinite column 
calculations, using equation ( 23 ) as a terminal boundary condi- 
tion. 

The analysis of chromatograph systems for pulse-type 
forcing functions has prompted consideration of the two ter- 
minal boundary conditions. The question arises as to how the 
use of a finite terminal boundary condition affects output 
prediction as compared to the Infinite column case when the 
system is forced by pulse- type functions. It Is desirable for 
the two predictions of column outlet behavior (z=l) to be sim- 
ilar so that the use of a priori estimates of Pe are valid In 
complicated models having the form of equation (22). 

In order to answer the above question and to establish 
the conditions under which a finite terminal boundary can be 
used to yield Inf inite column behavior at the column outlet 
(z=l), two relatively simple problems can be . considered: 


Case I: 


(l/Pe)(6V/Sz*) - = 6y/bQ 

(25) 

y(z,0) = 0; z>0 

(26) 

y(o.e) = <5(©); ©iO 

(27) 

11m y(z,e) = finite; ©>0 

(28) 

It 

o 

(29) 
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and j 

Case II ; j 

(l/l?e){<^y/dz^) - by/hz - = hy/bQ j ( 30 ) 

y(z,0) = 0; z>0 ■ j (31 ) 

y(o,e) = <((e): sio ( 32 ) 

6y(zQ,Q)/bz = 0; e>0, and arbitrary (33) 

Ra = 0 ( 3 ^) 


Case I considers the impulse response of the simple, 
one-dimensional, axial dispersion-convection model in an In- 
finite column. Case II considers the unit impulse response 
of the simple, one-dimensional, axial dispersion-convection 
model with the finite boundary condition. It is desirable to 
determine the conditions under which the tiio responses are 
equivalent. These conditions can be determined VJithout resort- 
ing to the comparislons of the analytical solutions for each 
case, through use of the method of moments. 

At a dimensionless length of unity, the column out- 
let, the Laplace transforms of the two solutions are^^*; 

Case I ; 

y(i,s) = exp[(Pe/2)-(arg)l (35) 

Case II ; 

y(l,s) = exp(Pe/2) ^[(Pe/2)+(arg)] exp[-(l-ZQ) (arg)] 

- [(Pe/2)-(arg)] exp^” (1 -Zq) (arg)]^^ 

^[(Pe/2)+(arg)] exp(Zo(arg)) 

- C(Pe/2)-(arg)] exp(-ZQ(arg) )^ 

* see Appendix E for details. 


( 36 ) 



where 


arg r= 


2 

(Pe/2) + Pe s 


(37) 


Each respective output curve can be characterized 
by Its moments. Two moments are considered here - the first 
moment about the origin and the second moment about the mean. 
The first moment about the origin gives the time of appearance 
of the mean of the output curve. The second moment about the 
mean gives the variance of the output curve. These moments, 
as has: been previously noted in Part 3» are directly obtain- 
able from the Laplace transform domain solution. The general 
relationships were given in equations (1^) through (I?)* Us- 
ing these relationships, the Case I and Case II transfer 
functions were analysed to yield: 

Mij = 1 (38) 

_ ^2i = 2/Pe (39) 

and 

= 1 + [exp(-Pe Zq) - exp(Pe - Pe Zq)] /Pe (40) 

Aajj = 2/Pe + exp(Pe - Pe z^) [4/Pe - ^z^/Pe ‘ 

- 4/Pe^] + exp(-2 Zq Pe)/Pe^ 

- exp(2 Pe - 2 Zq Pe)/Pe^ (41) 

If one considers the limit of the Case II moments as z^ becomes 
very large, the two results are equivalent, or: 

11m 1 

lim Mitt “ 2/Pe 


and 
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Table 3 suinmarlzes the results of parametric stud- 
ies of the two moments considered for each case. The errors 
in Case II versus Case I moments for 2^=1 are significant for 
low Peclet number* The error diminishes with increasing Peclet 
number. This confirms the qualitative conclusions of Friedly 
(1972) for high values of Pe. Table 3 also gives the value 
of Zq which, when used in Case II, will yield output character- 
istics the same as Case I output characteristics. This means 
that for a given Peclet number, application of equation (33) 
at the noted z^, will yield output characteristics at z=l that 
are, for all intents and purposes, the same as those predicted 

V 

by Case I. 

Table 4 presents some typical values of the Peclet 
parameter for several systems. For chromatographic systems, 
the range of the Peclet number is on the order of 5»000 to 
10,000. Thus in this research, it appears that use of the 
zero-derivative condition (equation (33)) at the column exit 
will not, cause serious problems. 

In conclusion, the comparlsion of the mean and var- 
iance for impulse responses at z=l for the two different bound- 
ary conditions has yielded guidelines which are useful when 
approximating infinite column behavior using a finite terminal 
boundary condition. The use of the criteria for general 
pulse-type forcing functions would yield results wherein the 
absolute errors between the two cases would be the same but 
the relative errors between cases would decrease. The guide- 
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TABLE 3 

Case I and Case II Comparislon Results 


Errors at 1.0 


Pe 

Absolute 

Error 

Relative 
Error , % 

' i /^ 1 j; ^ )UO 0 

Absolute 

Error 

Relative 
„ Error, ^ 

2 

0.4323 

43.23 

1.245 

124.5 

4 

0.2454 

24. 54 

0.3125 

62.9 

8 

0.1250 

12.50 

0. 0781 3 

31.2 

16 

0,06250 

6.250 

0.01953 

15.6 

32 

0,03125 

3.125 

0.004883 

7.91 

64 

0.01 563 

1.563 

0.001221 

3.91 

128 

0,00781 

0.781 

0.0003052 

1.99 

256 

0.00391 

0.391 

0.0000763 

0.976 

512 

0,00195 

0.195 

0.0000191 

0.489 

1024 

0.00098 

0.096 

0.000004? 

0.241 

2048 

0.00049 

0.049 

0.0000012 

0.123 

4096 

0.00024 

0.024 

0,0000003 

0.0615 

8192 

0.00012 

0.012 

0.0000000 

0.0 


Case I characteristics = Case II characteristics 





Safe Zq 

Safe Zq 


Pe 


[M\x “ n )* 

{U2X- )*«■ 


2 


9.791 

11 . 768 


4 


5.254 

6.021 


8. 


3.043 ■ 

3.328 


16 


1.978 

2.073 


32 


1.467 

1.490 


64 


1.223 

1.222 


128 


1.106 

1.044 


256 


1.050 

1.002 


512 


1.024 

1.0005 


1024 


1.011 

1.0001 


2048 


1.005 

1.00002 


4096 


1.002 

1.000005 


8192 


1.001 

1 . 000001 

* (-«.I 

- ) 

0 

1 

00 



** (Aji 

- ) 

IN 

»-* 

0 

I 

CD 
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TABLE 4 

Peolet Numbers for Four Typical Systems 


System 

Pe 

Reference 

Micro Gas Chromatograph Column 
(Water In Helium) 

233 

(Wilhite, 1966^ 

Typical Gas Chromatograph Column 
(Water in Helium) 

5622 

(Keba and Woodrow, 

1972) 

Typical Gas Dehydrator 
(Water in Helium) 

1777 

(Lashmet, 1973) 

Small Experimental Reactor 

155 

(Smith, 1970) 


(SO2 In Air) 
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lines developed here should also prove useful for models in- 
cluding other transport, mechanisms 0). When applicable 

to more complicated models, the method of analysis used here 
will give more definite guidelines for each specific situation. 



PART 5 

ORTHOGONAL COLLOCATION AS A NUMERICAL TECHNIQUE 

A* Motivation for Study of Orthopconal Collocation 

The complexity of the model formulation in Part 2 
necessitates the application of numerical approximation tech- 
niques to effect solution of the system of partial differen- 
^ tlal equations. A preliminary study of the widely prevalent 
technique known as Finite Differences has been made to ascer- 
tain whether or not this computational technique would prove 
suitable and effective for solution of the type of equations 
encountered in chromatograph system modeling. 

Finite difference approximations have predominantly 
been used in the analysis of partial differential equations. 

To obtain numerical solutions to partial differential equa- 
tions, one replaces the continuous variables with discrete 
variables. The relations between these discrete variables in 
the method of finite differences are called finite difference 
equations. The relationships are based on Taylor series rep- 
resentations of the dependent variable. The domains of the 
independent variables that are discretized form a system of 
grid points. Figure 5 shows a grid representation for the 
. transient analysis of a system with one spatial independent 
variable. The spatial dimension, z, is shown as being bounded 
and the time variable, 9, is shovm with no particular bound. 
The grid is a fixed grid; i.e., spatial discretizations and 
time discretizations are uniform for each domain. Note that 


3 ^ 




FIGURE 6 
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the value of z, the continuous space dimension is given by; 

2 == 1- { A Z ) j 

I 

where 1 refers to a particular spatial grid liolnt and Az is 

i ' 

the spacing between spatial grid points. Similarly, the 
value of 0, the continuous time variable is given by; 

e = j*{ ) 

where J refers to a particular time grid point and A6 is the 
interval between time grid points. 

For parabolic problems (as is the case for the sec- 
ond-order chromatograph system models), the two-level implicit 
method knovm as the Crank-Nicolson method is probably most 
popular and is well documented (Lapidus, 1962). In this method, 
the follov/ing approximations are made for the first and second 
spatial derivatives and the first time derivative: 


/s " J j+i 

(dy/6z;i — — ; + ; — ; — — 

] 2 (Az) 2 (Az) 


^ (Az)^ 


(AZ) 


(dy/de)ij= (yi.j+i - yij)/(*e) 


where the i subscript denotes a coordinate in the spatial 
domain and the J subscript denotes a coordinate in the time 
domain. 

Preliminary studies have been made applying the 
Crank-Nlcolson method to the problem: 



37 


(l/Pe) (^?y/5z^) - by/dz = by/dQ 

y(z,0) = 0; z-0 
y(o,e) = 0(e); e>o 
6y(zo#6)/i0 = 0; 0>O 

1 

' l 

Simulations were made with following conditions: 

1. 0(e) was a triangular-type pulse of duration 
0.01 and with unit area. This is quite a sharp 
pulse as far as typical chromatograph input 
pulses are concerned, hut it was used mainly in 
the interest of saving computer time. 

2. The Peclet was fixed at 8,000. 

3. The time increment, A0, was held at 0.0004 

4. The response was studied at z=0.05. This is a 
drastic reduction in the normal spatial coor- 
dinate studied, but, again, this was in the in- 

' terest of conserving computer time. 

5. The terminal boundary condition was applied at 

Zq = 0.20. 

6. The spatial Increment, Az, was varied in the 
following sequence: 

0.0002, 0.0004, 0.0010, 0.0025 
For spatial increment values of 0.0010 and less, the simula- 
tions were stable. However, when az was Increased to 0.0025. 
instability in the form of oscillation in the response was 
exhibited. The very small az required is directly attributable 
to the Pe value used. This instable az value is not quite as 
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small as the value that is predicted by the stability of 
Price, ( 1966 ) 

The simulations for spatial Increments of 0.0002, 
0.0004, 0.0010 gave reasonable results when compared to results 
convolving 0(6) with the analytical impulse response. The 
discrepancy between the analytic and numerical computations 
appeared in the magnitudes of each response point - the numer- 
ical results viere on the order of 20^ too low. This in turn 
affected the areas beneath the response curve for the numer- 
ical results - all areas were on the order of 0.80 as compared 
with the correct area of 1.0^ The area under the analytical 
response curve was 0,96 which is tolerable considering the 
sharp input. This discrepancy in response area can be resolved 
by adding additional parameters to the difference equations 
to yield an exact conservative relationship (Rogers, 1973) s 

[ System Input - System Output over*] 
the interval j to J +1 

/si i^\ 

where N is the total number of spatial points* This analysis 
was not performed because it was felt that the method already 
suffered from a more alarming feature - the high degree of 
spatial discretization v/hlch is necessary for the large Pe 
values encountered in chromatographic systems analysis. Ex- 
trapolation of the computing time required for the simulations 
performed yields an estimate of one to two hours of computer 
time (IBM 360 / 50 , F0RTRAN G) required for complete simulations 
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over the space interval (0, l.C"**)* The time would naturally 
Increase when broader input pulses are used. Similar conclu- 
sions on the use of finite difference schemes were reported 
earlier (Pfeiffer, I972), 

Because of the high degree of spatial discretization 
required by the finite difference method and the subsequent 
high cost of computer simulations, it was felt that further 
pursuit of finite difference formulations for problems sim- 
ilar to the above v;as not warranted in this investigation and 
that other techniques should be studied to determine if they 
would be computationally more expedient and desirable. 

Theory and Background of Orthop:onal Collocation 

A recent text (Pinlayson, 19?2) has dealt with sev- 
eral approximation techniques for the solution of the differ- 
ential equations which arise in the analysis of transport 
phenomena. A group of approximation techniques has been des- 
ignated the Method of Weighted Residuals (MWR). A subclass 
of MV/R is the Method of Orthogonal Collocation, This method 
has been successfully applied to several problems in the realm 
of chemical reaction engineering. Investigators In this area 
Include Ferguson and Pinlayson (1970), Pinlayson (197I), 
Vllladsen and Stevzart (I96?), Villadsen and Sorensen (I969), 
and Villadsen (197O), The purpose of this section is to pre- 
sent a summary of the theory behind the method. Discussion of 
investigations that concern general computational aspects will 
follow In the following section, . 
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The Method of Weighted Residuals approach to the 
solution of partial differential equations starts with a rep- 
resentation of of the dependent variable, y, by a finite sura 
of trial functions P 4 . An example might bei 

y(z, 0 ) = 0 Q(z.e) + ^l(^) (^ 2 ) 

4 . — » 

/a I 

where 0 q{z, 6 ) is a function which may be chosen to satisfy 
one or more boundary conditions. The functions Pj^(z) are nor- 
mally specified and the time-varying coefficients, a^(e), are^ 
determined in a manner to give the ”best*‘ solution of the dif- 
ferential equation. 

The next step in the MWR is to manipulate the differ 
ential equation such that one side, say the right hand side, 
of the equation is zero. Then,, the trial function expansion 
is substituted into the left hand side. This substitution of 
the trial function expansion into the manipulated differential 
equation forms what is termed the residual. Res. If the trial 
function were exact, the residual would be zero. In MKR, the 
coefficients, ai(9), are determined by specifying Tvelghted 
integrals of the residual to be zero; l.e.. 


Wj (Res) dV = 0; j=1.2, ... N 


(43) 


The choice of weighting functions, Wj, determines 
what class of HV/R is to be applied. For the general colloca- 
tion method, the weighting functions are chosen as displaced 
Pirac delta functions: 

Wj = <i(z-Zj);j=l,2, ...N ( 44 ) 


N 



Substitution of equation (^4) Into equation (43) gives the 

j 

} 

result of forcing the residual to be zero at N specified col- 
location points* As the degree of approximation is increased, 

i 

the residual will forced to be zero at an increasing number 
of points in the spatial domain and the trial function should 
converge to the true solution within a given accuracy. 

Within the class of collocation methods is the sub- 
class of orthogonal collocation. The distinguishing feature 
of this method is that the trial functions* Pj^(z), are chosen 

as orthogonal polynomials defined by the following relationship; 
b 

J*w{z) Pi(z) Pj(z) (45) 

cZ 

where ja, is the interval of orthogonality, w(z) is a pos- 
itive weighting function on [n,b^ , is a scale factor, and 
is the Kronecker delta. The group of polynomials defined 
by equation (44) is said to be orthogonal on the interval 
with respect to the weighting function w(z). 

The N collocation points are chosen as roots to 
Pjj(z), which is the polynomial of the next highest order in 
the trial function expansion, the highest being P|j^i in equa- 
tion (42). The basis for choosing the roots of the polynomial 
as the collocation points instead of equidistant points in the 
interval of interest can be found in the theory of polynomial 
interpolation. Several results, as documented by Lanczos (1956) 
are summarized here; 

1, Polynomial expansions are justified due to the 



fundamental theorem proved by Welerstrass in 
1885 '^hlch establishes that any continuous 
function In a finite Interval can always be 
approximated to any degree of accuracy by 
finite power series. 

2. The V/eiers trass theorem does not imply that an 
approximating polynomial can be obtained by us- 
ing equidistant points. This behavior was stud- 
ied by Runge In I 9 OI who showed that equidistant 
interpolation of some very simple analytical 
functions could In certain regions yield very 
erroneous results which did not disappear v:ith 
increased points. This behavior is termed the 
”Runge phenomenon.” 

3. The difficulties which occur with equidistant 
Interpolation disappear vfhen the zeros of the 
first neglected polynomial in the polynomial 
approximation are used as the interpolation 
points. Hov?everp this Introduces the need to 
loiow the roots of the particular polynomial. 

C. General Computational Aspects of Orthogonal Collocation 

The solution of parabolic partial differential 
equations using orthogonal collocation requires several steps 
which are Independent of the particular equation under con- 
sideration, This section presents two formulations which are 
theoretically equi1r6-lent. but which differ in computational 
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and coding advantages. The first formulation, although some- 
what more complex from a coding point of view, will be shovm 
to be superior for computations. 

A trial function has been proposed, Finlayson (1972, 
p, 105 ), for second order systems on the spatial Interval [o,l] • 
For transient analysis, the trial function Is of the form: 

N 

^ y(z,e) = f(e) + s(e) z + z(i^z)^ ai(o) p±^i(z) (^6) 

i==i 

The above equation has N+2 unknowns: The functions f(9), g(6), 

and aj^(9) p 1=1, N, These are determined by the boundary con- 
ditions at z=0 and z=l and by performing collocation at the K 
roots of Thus one has a set of N+2 points: 


zi = 0 

^N+2= ^ 


and N+1 ; the roots of 

Now, if one were to construct the approximate solu- 
tlon at these N+2 points, a matrix equation would result: 

y(zi,e) 1 Zi z^(1--Zj^)Pq(z^) ...Z3^(1-Z2^)P^(zjl) f(9) 

y(z2*e) 1 ^2 Z2^^“22)^o(^25 • • 

• » • 4 ^ 


y(zi,e) 


y(z2*^) 

• 


• 

• 

y(zH+i-®) 


y(zN+2-®) 



1 %+i 


[y ( ® I ^ ^N+ 2 ^ ^0 ( • 

Now define the following quantities: 


^ ^0 ( ^ ajj ( G ) 




A 

y «s 


R 


f = 


f y(zi,e) 

yCzg.e) 


y(zN+i.0) 

y{zn+2*®^ 


1 zi zi(1-zi)Pq(zi) 

1 zg Z2(l-Z2)Po(z2) 

• • 

• • 

i zn+1 

f(e) 
s(e) 
a^Ce) 

ajj(e) 


Zl(l-zl)P„-|(zi) 

Z2(l— Z2)^j-i( 22 ) 


(48) 


(49) 


(50) 


Use of equations (48) > (49) and (5f) reduces equation i4?) to 


the more compact form: 

Z « i f (51) 

The spatial derivatives may he expressed in a similar form: 


6y/6z = U f 


(52) 


where 


^y/6z^ = H f 


^y(zi .e)/^2 

^y(Z;i,9)/dZ 


Jyy/6z = 


hy{z^^,9)/tz 

by(z^,j,o)/6z 


(53) 


(54) 
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and 




N j.i 
[^j J.2 
M j.i 


iy{zi ,e)/hz^ 
^y{zz^e)/bz^ 


Mz^^^e)/bz^ 


^ 0 ; j=l,N+2 
= 1; j=l,N+2 

« 2j(l-Zj) dP3^^^(Zj)/d2 

+ (1.2zj) Pi.3(z_j) ; J=1 


N+2 


1 « 3 p N+2 


H i.i 
N j.i 


= [ aa ] >2, - Os j = i . N+2 

= Zj(l-Zj) d^Pj^_3(Zj)/dz^ 

+ 2(l-2zp dP^^^(zj)/dz 


- 2Pi_3(Zj) 


; ;J=lpN+2 
i=^3*N+2 


(55) 


(56) 


(57) 


The time-varying vector f may be eliminated from equations (52) 
and (53) by premultlpling equation (5^) by the inverse of R» 


R“ 


or? 


and 


-1 

f : = R y 


6y/bz = gl y 

dV/dz^ = R2 y 


(58) 

(59) 


Equations (58) and (59) thus yield expressions for the first 



and second spatial derivatives at the N+2 points in terms of 
the solution at the K+2 points. 

Alternative to the formulation of above Is a form- 
ulation which Is presented by Flnlayson (197^* PP* 105-106). 

Expansion of equation (46) yields an (N+1) order polynomial: 

N+1 

y(z,G) = f(e) + S d4(e) (6o) 

1=1 


Writing the approximate solution at the N+2 points yields a 
matrix equation similar to equation ( 51 )s 


y « Q d 


( 61 ) 


where: 


Q ^ 


21 


2 

2l 


.^N+l^K-fl** • 
2 


N+1 


N+1 

'N+1 

K+1 

'N+2 


( 62 ) 


d 4 


f(e) 

di(G) 

d2(©) 


d (e) 


( 63 ) 


The first and second spatial derivative vectors can be writ- 
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where : 


i>y/dz 

= Ql d 


« Q2 d 

W « 

»5 (i-i) 

Nji 

= (i-l)(l-2) z 


j=l,N+2 
1=1, N+2 

J=l.N+2 
'* 1=1, N+2 


(64) 

(65) 


( 66 ) 

( 6 ?) 


As in the first formulation, the time-varying vector, d, may 

be eliminated from equations (64) and (65) by pre-multiplylng 

1 

equation (61) by the inverse of Q, 5T » 


d 0“^ 2. 

and 

dV^z^ = S§ z 


( 68 ) 

(69) 


Thus, equations (68) and (69) give expressions which are Idrn- 
tlcal to equations (58) (59)- The matrix product. Rl 

is equivalent to and R2 is equivalent to g2 • 

Since the computations of g, and §2 only require know- 

ledge of the collocation points and not knovfledge of the par- 
ticular polynomial coefficients being considered, one might 
conceivably prefer the second formulation. Both formulations 
require the computation of the Inverse of an (N+2) square 


matrix. 



.1 


kQ 

Computationally^ it Is desirable for the matrix 
being Inverted to be well-conditioned with respect to inver- 
sion. An analysis was made comparing the inversion qualities 
of the matrices R and The ease of Inversion is measured 
by the condition number of ^ and g respectively and with the 
number of decimal digits which are left unchanged following 
iterative improvement of the initial Gauss-Jordan reduction 
of each matrix, Stewart (1973) discusses the problem of Ill- 
conditioning and the use of iterative improvement in matrix 
Inversion. Table 5 compares the inversion characteristics of 
R and ^ for increasing N. The condition numbers cited are 
lower bounds on the true condition numbers relative to the 
% norm*. Appendix F shows how the lower bound and upper 
bound on the condition number Is computed. Except in the 
analysis of for (N+2) - 22, there were no practical dif- 
ferences in the lower and upper bounds. 

Table 5 indicates that the matrix g is well-conditioned 
with respect to inversion using the double precision word 
length available on. the IBM 3^0/50 computer. In all cases, 
the computation of the product £ yielded a matrix whose 
off-diagonal elements were less than or equal to 10“* • The 

table also shovs the progressively poorer conditioning of Q 
with respect to inversion* The (26 x 26) case is so Hi- 
tt the hi norm of an (n x n) matrix A is defined as: 


Ll norm(A) 


s= max 






n 



k? 


TABLE 5 

( 

Comparlslon of Conditioning of E 
Matrices with Respect to Inversion ! 


Matrix 1 

Size 

Lower Bound 
on Condition of R 

IDGTg 

Lower Bound 
on Condition of 

IDGT^ 
Q ^ 

{ 3 

X 

3 

T“ 

0.120 

X 

100 

i5 

0.240 

X 

100 


( 

X 

4 

) 

0.328 

X 

10^ 

15 

0.149 

X 

103 

15 

( 5 

X 

5 

) 

0.739 

X 

10^ 

15 

0.944 

X 

10^ 

15 

( 6 

X 

6 

) 

0.142 

X 

10^ 

15 

0.591 

X 

lo'^ 

14 

( 7 

X 

7 

) 

0.243 

X 

10^ 

15 

0.366 

X 

10 ^ 

14 

( 8 

X 

.8 

) 

0.384 

X 

10^ 

15 

0.225 

X 

xo^ 

13 

( 9 

X 

9 

) 

0.571 

X 

10^ 

15 

0.138 

X 

lo7 

13 

(10 

X 

10) 

0.812 

X 

10^ 

15 

0.840 

X 

10^ 

11 

(11 

X 

11) 

0.111 

X 

lo'' 

15 , 

0.510 

X 

10® 

11 

(12 

X 

12) 

0.148 

X 

10^ 

15 

0.309 

X 

109 

11 

(lij- 

X 

Ik) 

0.244 

X 

10^ 

15 

0.112 

X 

10^^ 

9 

(18 

X 

18) 

0.545 

X 

10^ 

15 

0.145 

X 

10^^ 

6 

(22 

X 

22) 

0.103 

X 

10^ 

15 

0.177 

X 

10^^ 

2 

(26 

X 

26 ) 

0.179 

X 

10^ 

15 

0.907 

X 

10^® 

0*^ 


* IDGT Is the approximate number of digits in the Inverse which 
were left unchanged after iterative Improvement 
There was no convergence in the iterative improvement. The 
upper bound on the condition of g was 0.202 x 10^^ based on 
the ”best‘* . 

Note t Subscripts R and Q on IDGT refer to inversion of R and 
Q respectively. ~ 



50 


conditioned that inversion using the available computer is 
computationally Impossible. Based on condition numbers and 

accuracy, either formulation is acceptable fob (N+2) ^ 5; 

I 

while for (N+2) > 6 , the first formulation is preferable « 

Ife should be noted that Plnlayson(1972, p»35) 
alludes to this problem but the comment is quite qualita- 
tive and somewhat obscure: 

“The orthogoiiallty of the polynomials 
gives computational advantages, although the 
same approximation can be expressed in terms 
of powers of x, if the computations can be 
done accurately enough^. 

The preceding analysis used the roots of the so-called 
shifted Legendre polynomials. These are defined by equation 
(45) if one lets a=0, b=l , and w( 2 )=lo The polynomial coef- 
ficients were computed using the relationships of Villadsen 
(1970). Figure 6 shows the behavior of the first four of 
these polynomials. The roots were computed by ' shifting the 
abscissas from Gaussian quadrature formulae, available in 
Abramowitz and Segun ( 1965)0 Love (I 966 ), and Stroud and 
Secrest (I 966 )* 

Although most of the problems solved by others us- 
ing orthogonal collocation have not required over 12 collocation 
points, the results of this section point out a computation 
disadvantage of the second formulation which appears at a 
fairly small degree of discretization and gets progressively 
worse* The first formulation requires some additional Infor- 
mation but successfully circumvents the problems Inherent in 




the second formulation. Of course, the precision capabilities 
of the computer used for computations must be taken into 
account also. 

In addition, these remarks carry over into problems 
where it is chosen to use polynomials in the squared spatial 
dimension. For example, a trial function which has been pro- 
posed for solution of a catalyst modeling problem isi 

N 

y(yi,0) ~ h(9) + (l-A®') ^ a^(0) ( 70 ) 

where h(0) is determined by the boundary condition at A=1 
and the boundary condition ^y(0,©)/^ ^^0 is satisfied by the 
trial function. Subsequent expansion and formulation at the 
respective collocation points yields a matrix to be inverted 
and the procedure of retaining polynomials within a coefficient 
matrix similar to R is favored over a formulation yielding a 
matrix similar to Q for the reasons previously listed. 



PART 6 


APPLICATION OP ORTHOGONAL COLLOCATION TO A 
TRANSIENT, LIFFU SI ON-CONVECTION MASS TRANSFER 

PROBLEM 

The use of orthogonal collocation as a technique 
for solution of pulsed, distributed systems, the chromatograph 
system being only one such system, is an area which has not been 
documented in current literature. Because of this lack: of c con- 
tribution in this area, guidelines for effective use of this 
method must be established and documented. 

In this section the general aspects of orthogonal 
collocation enumerated upon in the previous section will be 
applied to solve a simple, transient diffusion-convection mass 


transfer problem; 

(1/Pe ) (^y/6z^) - by/dz = by/bO (71) 

. y(z»0) = Oj 2>0 (72) 

y(o,e) « 0(e)t QzQ (73) 


Motivation for the study of this problem is several 
fold. First of all, the problem has a direct analj^tlc solution, 
therefore giving a result, useful for comparlsion. Secondly, 
the problem possesses characteristics of more complex models. 
Thirdly, successful application of orthogonal collocation 
should give guidelines for subsequent applications. 
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The preceding analysis has "been conducted based on 
a spatial interval of [o,l] as the interval of orthogonality 
for the orthogonal polynomials In the trial function expansion. 


However* as v^as shown in Part application of boundary con- 
dition (?4) should be made at some point 2 ^ which should be 
different than unity depending on the value of the Peclet 
, number. To avoid derivation of additional polynomials orthog- 
onal on the interval [0,Zo] determination of the required 

roots, the above problem may be rescaled in the spatial domain 
by the following change in variable: 

~ ( t ) z 

Therefor© 

(l/dz) s= (1 /Zq) (l/dzj^gy) (75) 

a/hzf = (1/Zo)^(1AznEw)^ (76) 

Use of equations (75) and (7^) and deletion of the subscript 
“NEW*’ yields the rescaled problem: 


(l/Pe)(l/Zj,)^(5*y/6z^) • 
y(z,0) = 0; Z>0 
y(0,e) = 0(0); 050 
6y(l,0)/dz = 0; 0>O 


(l/z^)by/6z «= hy/bQ (77) 

(78) 

(79) 

(80) 


Where one was concerned with the dimensionless length of uni- 
ty in the old coordinate system, one Is now concerned with 
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the dimensionless length of (1 /zq) which now corresponds to 
the outlet of the bed. 

As a variant of the favored formulation of Part 5» 
one may represent the solution at the N collocation points, 

[22, ZN+l] I asj 

y(z2 «®) 1 Z2 

= . f(e) + 

y(2jg+l,e)J [ij [^K+l 

Z2(1-Z2>Po(^2> 

# 

« 

Formulation of the problem in this manner reduces the size of 
the matrix which must be inverted from (N+2)x(N+2) to NxI'J, al-» 
though Increasing the coding effort. Application of equation 
(81) to the above problem reduces the distributed system to a 
set of N ordinary differential equations represented by: 

i « w i 0 (e)); £(0) = 0 (82) 

Manipulative details and full matrix definitions for this 
problem are provided in Appondices G and H. It should be noted 
that the function, w(z), which appears outside the summation 
sign in the general trial function (Appendix G) has been taken 
as vr(z) = z(l-z) as previously seen in equation (46). This 
specific form is used exclusively in this investigation for 



g(e) 

. .. Z2(1-Z2)%^i(z2^ 

* • • • • ( 81 ) 

... ajj(0) 
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all of the interparticle approximations. However, one might 

I 

possibly specify the form of w(z) relative to the types of 

polynomials used. That Is, for a given function w(z), one 

i 

might specify the polynomials such that the columns of the 

j 

NxN coefficient matrix In equation (81) become orthogonal. 

Thus, one would be taking advantage of the orthogonality 
properties of the specific polynomials rather than using an 
arbitrary polynomial set. For example, If one used w(z) = 
z(l-z) In the trial function (as is shown In equation (81)), 
the columns of the NxN matrix could be made orthogonal by 
defining the polynomial set byi 

|z^(l-z)^Pj^(z) Pj(z) 

The Inner products of the matrix columns would represent the 
discrete form of this Integral. 

The merits of utilizing the orthogonality properties 
of the specified function set has not been established. As 
will pointed out in Part 9, Discussion, the undertaking of such 
a study requires computing capabilities (precision) to deter- 
mine roots of polynomials which may be ” uncommon^* and not 
tabulated to a large number of significant figures. For these 
reasons, the choice of polynomials In this investigation was 
dictated by the availability of the high precision roots. 

An elgenanalysis* of the matrix W In equation (82) 

*A computer "^"ogram listing is given as part of Appendix H. 

This program performed all the manipulations and computations 
documented In Appendices G and H as well as the elgenanalysls . 
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WQS made for Peclet numbers of 1, 10, 100, 1000, and 10000. 

I 

The number of collocation points, N, Kas varl<^d In the se- 
quence 3, 4, 7, 15, and 21. The shifted Legeridre polynomials 
and roots were used In the analysis. The valqe of Zq was 
held at 2.0 for all cases. 

For all cases, the eigenanalysis yielded eigenvalues 
with negative real parts Indicating a stable set of ordinary 
differential equations. This result contradicts the results 
obtained earlier (Woodrow, 1973) for Pe=l. The difference 
between the analysis lies in the sequence of manipulations 
and computations made In arriving at the matrix. In this paper 
W, to be analysed. Although the approaches are equivalent 
theoretically, the comput potions produced different. The approach 
detailed in Appendix G is therefore favored. 

While stability is Indicated by the negative real 
parts of the. eigen vaiUp s', an oscillatory behavior was indi- 
cated by the presence of Imaginary parts for a majority of 
the eigenvalues in each case. The magnitudes of the imagi- 
nary parts Increased with Increasing Pe. Therefore, it was 
reasonable to expect that simulations using the orthogonal 
collocation technique would exhibit some degree of oscilla- 
tion depending on how the modes of the matrix, W, were 
coupled. 

Various simulations have been performed for this 
problem. Table 6 sumarizes the different cases considered in 
this investigation. The method by which the set of differential 
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TABLE 6 


Summary of Orthogonal Collocation Slmulat^ions 


for the Trans 

lent, Diffusion-Convection 

problem 

0(e) =/unit 
Pe 

rectangulai' pulse of five 
N ' Figure Reference 

times duratloi^ 
Execution 

32*0 

3 

Figure 8 

( sec. ) 

23.55 

32.0 

7 

Figure 9 

102.05 

32.0 

15 

Figure 10 

ll 6 i )-.32 

10000.0 

15 

Figure 12 

301 , 62 

10000.0 

21 

Figure 13 

891.10 

0(e) = Actual input 

data. Figure 14^ 


Pe 

N 

Figure Reference 

Execution tlrae«» 

10000.0 

3 

Figure 16 

{sec. ) 

. 23.90 

10000.0 

7 

Figure 17 

57-10 

10000.0 

15 

Figure 18 

419.25 


^ Double precision computations using F0RTRAN G on IBM 360^5^* 
Integrat ons terminated at ten time units. 

«^*«-Double precision computations using F0RTRAN G on IBM 380/50; 
integrations terminated at twenty time units. 

Note t For all cases = 2.0; responses for all collocation 
points outputed at each time Increment. 
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equations was solved Is based on an extrapolation treatment 
(Bullrsch and Stoer, I966). This algorithm appears in a sub- 
routine (DREBS) which Is part of the IMSL scientific subrou- 
tine library (IMSL, 1973» P- LREBS) currently available on 
Rensselaer’s IBM 360/50 computing facility. The computations 
were made entirely in double precision using the F0RTRAN G 
compiler. Table 6 shov/s how execution time for the simula- 
tions was affected by N, the number of collocation points. 

The method used to integrate the equations is not 
too well kno™. In their paper, Bulirsch and Stoer showed 
the superiority of the extrapolation treatment over a Runge- 
Kutta and Adams-Koulton-Bashforth methods, A compariaion v?as 
mads between the subroutine DREBS, the IBM SSP (IBM, I968) 
subroutines for a fourth order Runge-Kutta and Hamming Predic- 
tor-Corrector method, and a simple Euler Method, The com- 
parision was based on the N=3 simulation for this system. For 
the same error criteria, it was found that the Euler method 
was significantly slower than the extrapolation treatment and 
while the Runge-Kutta and Predictor-Corrector methods used a 
larger step size than the Euler step size, the step was still 
much smaller than the extrapolation treatment and hence was 
computationally slower. This result agrees with Bulirsch and 
Stoer for the problems that they considered. 

Although the eigenanalysls indicated that the system 
of ordinary differential equations was stable, a closer exam- 
ination (made near the conclusion of this investigation) of the 
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computed eigenvalues Indicated that the system became increas- 
ingly •’stiff” with Increasing approximation order. Table 7 
shows this behavior and the behavior as affected by Feclet 
number. The parenthesised ratios are the absolute value of 
the largest eigenvalue real part to the absolute value of the 
smallest eigenvalue real part. The result of the indicated 
division is shown beneath each ratio. Using these ratios as 

I * • V ■ 

measures of stiffness i one can see that the eq.uation set is 
f stiff er” for low values of Peclet number and hence a smaller 
Integration step is required when the large eigenvalue response 
mode contributes to the solution* This would have the ^effect 
of increasing computation time with low Peclet niunber for a 
given order of discretization. This deduction is consistent 
with the increased computation times for the low Peclet number 
versus the high Peclet niimber simulations for the unit rec- 
tangular forcing pulse (see Table 6). Although the Bulirsch 
and Stoer extrapolation technique may not be particularly v;ell 
suited for the “stiff” system* it \:as used as the method of 
integration. In retrospect, another algorithm might have been 
better - perhaps a variable order Adams method (Hull et.al. , 
1972 ). V/ithin tho IMSL library* the subroutine DVG0ER (IMSL, 
1973 # P.DVG0ER) using Gear’s implementation (Gear, 1971a, 1971 1>) 
would bo a likely candidate for use. This situation could 
form an additional area of analysis - whether orthogonal col- 
location approximations produce, in general, stiff ordinary 
differential equation sets and what integration algorithm can 



TABLE 7 


Largest to Smallest Eigenvalue Ratios for 
the Orthogonal Collocation Discretization 
of the Simple, Dlf fusion-Convectlon Problem 


N. 



Pe 




1 

10 

100 

1000 

1 0000 

3 - 

/22.26^ 
' X . Zii ' 

= 17.4 

(liM) 

^1.19 
= 2.92 

= 1.63 

/ 1 . 46 \ 
'0.90 

= 1.63 

'o 7 B 9 
= 1.63 

7 

- 

. = 2 ? 4 . 

/ 44 . 1 \ 

'2,52' 
= 17.6 

( 10 .- 45 ) 
' 0.95 
= 11.1 

(2^) 

'0.57' 

= 4.53 

'o. 53 \ 
= 6,2 

15 

= 4700.0 

(6M) 
^ 2.70' 

= 236.0 

' 2,72' 

= 36.5 

= 40.3 

= 127.1 

21 

, 21 892 . 0 ^ 

< 1:25^ 

= 17100.0 

.2272.OV 
' 2.70' 

= 840.0 

= 82.0 

( 61 . 5 -) 
u.oo ' 

= 81.5 

(40 ,. 24 ) 
^ 0,30' 

= 136.0 


Ratios are ( max |^real | / ^^^1 ^^eall ^ 
where Vs are the eigenvalues of W In 



bd used most effectively for simulation purposes. 

Eoturnlng to Table 6, the first set ;of cases involves 
the use of a unit rectangular pulse of five time units dura- 
tion as the forcing function, 0(9). Within this first set 
of cases, the Peclet number was set at 32. The exact response, 
computed by convolution of the impulse response with the input 
is given in Figure 7. Figures 8, 9, and 10 show the orthogonal 
collocation approximations to the response for 7. and 15 
collocation points respectively. Note that all responses are 
for 2=0.5 and Zo=2,0. Hence, the responses correspond to the 
behavior at the bed outlet. This Is a convenient scaling of 
the problem because for the odd order approximations used and 
the shifted Legendre polynomials, the point 2=0,5 is always 
a root and hence collocation point. The response for 3 col- 

I 

location points shows several objectionable oscillations com- 

* ■ 

pared to the exact solution, although qualitative character- 
istics are well approximated. The response for seven collo- 
cation poi its exhibits several oscillations of much smaller 
amplitude and duration and the exact response is Increasingly 
well approximated. The response for fifteen collocation 
points is very close to the exact response and, within the 
accuracy of the plot. Is virtually Identical. However, the 
digital response did show small oscillations when the response 
“should have been” zero. 

Subsequent to the studies for Pe=32, it was decided 
to use a Peclet number more representative of the magnitude 



0 . 20 -- 

0.16- 


0 . 12 - 

y(O.5,0) 

0.0 8 - 


0.04 + 


0.0 



SIMPLE MODEL 
EXACT RESPONSE 
Pe = 32 
2=0.5 

2o=2,0 

T = 5.0 


8 .0 


-H 1— 

9.0 10.0 


-0.044- 



FI6URE 7 


I 






0.20 

0.16 

0.12 

y(O.5,0) 

0.08 

0.04 

0.0 

-0.04 



FIGURE 10 


COLLOCATION RESPONSE 
N = I5 

Pe = 32 
z = 0.5 

Zo“2.0 

T»5.0 


6.0 7.0 8.0 9.0 10.0 


Os 

OS 



67 


encountered In the chromatographic system, namely Pe« 10000* 
Figure 11 shows the exact response. With the large value of 
Pe, the character of the parabolic partial differential equa- 
tion becomes increasingly hyperbolic. The exact response is 
effectively the translated Input pulse vilth the corners slightly 
roimded and dispersed. The orthogonal collocation approxi- 
mations to the exact response sire given in Figures 1.2 and 13 
for N -15 and N=21 respectively.- Both approximations exhibit 
numerous moderate amplitude oscillations. This shows the dif- 
ficulty that the orthogonal collocation technique has in approx- 
imating functions with sharp, almost dlscontinsous behavior. 

For this situation, a high degree of discretization would be 
desirable. An attempt to' generate collocation matrices for a 
thirty-first order simulation was made.. This attempt was not 
successful because the matrix W showed instability in the form 
of positive eigenvalues. This result, which would completely 
reverse the trend of convergence to the solution with increas- 
ing N, was attributed to accumulated round-off errors in eval- 
uation of the coefficient matrix because the required preci- 
sion to carry the higher order polynomial coefficients becomes 
greater than the precision capability of the computer being 
used (IBM 360 / 50 ). 

Upon discovery of this weakness of the orthogonal 
collocation technique and the computational constraint of the 
IBM 360 / 50 , effort was directed to use of a “less” sharp 
forcing function in conjunction with the high Pe value. The 
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chosen forcing function was actual chromatograph Input data 
(normalized with respect to the dimensionless time, ©) shown 
In Figure 14, This Input data corresponds to an Injected 
Pentane sample at 150 C. studied previously (Keba and \7oodrow, 
1972 ). This forcing function is used for the second sot of 
cases listed In Table 6, The exact response, again computed 
by convolution of the impulse response with input is shown 
In Figure 15* Figures I 6 , 1?, and 18 show the orthogonal col- 
location approximations to the exact response for 3. 1 % and 
15 collocation points respectively. For purposes of numerical 
Integration, the input function was interpolated using cubic 
spline functions. Another IKSL subroutine, ICSICE (IMSL, 1973t 
p. ICSICE), was used to compute the interpolation coefficients. 
Again, all responses shown are for z=0.5 and Zq~ 2.0. The 
result for N=3, Figure 16 , exhibits an oscillatory behavior 
and gives a good qualitative representation of the true response. 
The result for N=7, Figure I 7 , exhibits a better approximation 
with reduced oscillations. The result for N=15* Figure 18, 
gives virtually the same result as the exact. Again, oscil- 
lations are still present in the digital results but are of 
very small magnitude. 

The preceding results offer several conclusions as 

to the usefulness of orthogonal collocation for the system 

» 

under consideration! 

1. Orthogonal collocation greatly reduces the degree 

of spatial discretization required for numerical 
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stability as compared to a finite difference dis- 
cretization, I 

2. High values of Peclet number combined with very 

j , 

• . j 

sharp forcing functions; rectangular pulses, 

require a degree of approximation which may become 
limited by computing capabilities* 

3# The use of smoother pulses in the cases of high 
Peclet number allows a very good result for N=15 
and a very reasonable result for N=7* 




PART 7 


APPLICATION OF ORTHOGONAL COLLOCATION TO A 
RATE OF ADSORPTION LIMITED IN TER- INTRAPARTICLE 

MODEL 

Subsequent to the studies of the previous section, 
attention was directed to applying orthogonal collocation to 
a model which was more complex and vjhlch might be used, under 
' certain conditions, as a viable model for a chromatographic 
system with porous packing material* If one considers the 
inter- Intraparticle model represented by equations (2), {3)*. 
(4), and (5) and considers the case where ^®A — 

the following model may be deduced (see Appendix A)j 

(l/Pe£)6^y/^z^ - hy/6z - Kj^uCy-y^) = ^y/^e (83) 

(l/Hj)6x^/de == Nj^u(y-y^) (84) 

y« == m x^ (85) 

For analytic solutions, the applicable Initial and boundary 
conditions are: 

y(z,0) = 0; z>0 (86) 

Xa(z.O) = 0; Z>0 (8?) 

y(o,e) = (i(e)j e^o (88) 

11m y(z,6) = finite! e>0 ( 89 ) 

2*^ CO 

Equations ( 83 ) through ( 89 ) form what Is termed the Rate of 
Adsorption Limited Inter^Intrapartlcle Model* For purposes 
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of numerical solution, the terminal boundary condition, 
equation (89)* is replaced by the finite terminal boundary 
condition previously applied* And the forcing function, 

(f(0), is replaced by a finite width pulse, 0(0)* 

This model is mathematically equivalent to an inter 
particle with finite rates of mass transfer which was pre- 
viously considered (Keba and Woodrow, 1972), By analogy, the 
Laplace transform of; the time domain response is (column 
outlet): 

y(l*s) = y(0,s) exp [(Peg/2) - arg] 

where , . . 

arg « \j Peg ^( s+NjjymEj ^ s+Nj^yiaRj )+Pejg/4 

+ Nj^y(l-mRj)J 

Prom equation (88), y(0,s) = 1* Substitution and inversion 
gives the Inlt Impulse response for the model, 

y(l,e) « o (y^ + y 2 > 

where: ' 

o «= (I) V^Pse/^) • exp (Peg/2 )» exp (-Nj^umRj-e) 

71 = dV^ ) exp [•(Pej2/40)-(Pe£;0/4)-N£u0 +N£umRxJ 

& • 

2 i \/NRUinRi(e-x)xJ 1 
72-2 NhuoHi I [2 ^ymRi ( 0-xIxJ vT 

• exp[-((Peg/4x)+(Pegx/4)+Ngy(l-mRj)x )] • dx 
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For inputs other than the unit Impulse, numerical convolution 
18 used to generate solutions. Using the techniques documented 
previously (Keba and Woodrow, 1972). a solution may be compu 
ted directly. The exact solution which Is subsequently pre- 

4- A this Dreviously documented technique* 

eented was computed using this pr 

This model has two parameters. mBi and Ngu. which 
are not estimable a £rlori and require determination via an 
appropriate curve fitting technique. In the example that 

th. 1= «.«.» t. th. 

.. th, m V.I., d,t.r.ln,a pr=vl.»»ly (Kh. -ooto". 

%u. 1- 

of th, output tot. "1th that "hlth 1. pro- 
toot. d W th, .od.l add,d to th, input tot. v.rl«.o, .to- 


tlon (20) )• 
section is 
102 column 


The data set that Is to be considered In this 
that for Acetone at 100 C. taken on the Chromasorb 
(Keba and Woodrow, 1972). The parameter mBi Is 


taken to be 0.029 and the parameter KHublsle.stlmated to be 
87.0. Figure 19 shows how %u «»s determined and for Qompa 
islon shows an equivalent relationship for the model devel- 


oped in Part 2. The plot shows that the neglecting of the 
diffusion (intrapartlcle) and mass transfer effects requires 
a smaller %o to give the same predicted output variance. 


Hence, the diffusive and mass transfer effects (primarily 
diffusive due to the high NtoG “lumped" Into 


the rate of adsorption parameter, Nbu* 
are the same as Indicated In Table 1. 


Other parameter values 
Figure 20 shows a plot 
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of the diraensionless time normalized Input data* Figure 21 
shows the exact response for this problem and the given param- 
eters computed using the analytic Impulse response and numer- 
ical convolution. 

Application of orthogonal, collocation to this model 
results in a set of 2N ordinary differential equations, where 
N Is again the number of collocation points. Appendix I gives 
the details of the orthogonal collocation approximation treat- 
ment of equations (83)» (84), and (85). The model is reduced 
to the following set of 2N ordinary differential equations: 

(90) 

(91) 

where 0(9) is again the forcing function and ^ “the 

(Nxl) vectors of the compositions and equilibrium compositions 
at the N collocation points, respectively. The matrices In 
equations (90) and (91) are fully defined In Appendix I. 

Several simulations using this orthogonal collocation 
approximation have been made. Table 8 gives a summary of 
these computations. The entries in this table are not in 
strict chronology. The following paragraphs fully present : 
these results and document the chronological details. In all 
oases, the system of ordinary differential equations was in- 
tegrated by using the previously referenced IMSL library sub- 
routine, DREBS. However, all simulations were made in single 


I = gi ^ + W2 jr* - W 10(0) } ^ 

(z “ 2*); r^co) = 0 
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TABLE 8 


Summary of Orthogonal Collocation Simulation 
Times for the Rato of Adsorption Limited,! 

Inter-Intrapartlcle Nodel l 


COMPUTER 
IBM 360/50^ 
IBM 360/50^ 
CDC 7600^ 


N = 3' 
(Figure 22) 


6*49 min. 


N = 7 

(Figure 23J 
25.05 min. 

23.53 min. 

0.316 min. 


N « 15 
(Figure 2^y 

110.22 min. 


N = 21 
(Figure 257 

390. min. 

(estimated) 

350. min. 

(estimated) 

2.87 min. 


All execution times are for single precision integration up to 
90 units of dimensionless time. 

^F0RTHAN H, output at every Integration step. 

^F0RTRAN H, output at approximately every 0.25 units of 
dimensionless time. 

^Output at approximately every 0.25 units of dimensionless 
time. 
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precision while some matrix addition and subtractions were 
made In double precision* A listing of the program used for 
these computations is given in Appendix I* As noted in 
Table 8, some responses v;ere computed on the IBM 3^0/50 
computing facility at RPI whereas others were made at the 
CDC 7^00 computing facility at Combustion Engineering* Inc** 
Windsor* Connecticut* As in Part 6, the spatial dimension 
was rescaled. The plotted results are for the collocation 
point 2=0. 5» vrlth 2^=2. 0* Again this corresponds to the bed 
outlet* 

Figure 22 shows the simulated response for N=3* 

This approximation Is extremely qualitative as compared to 
the exact solution. The only correct prediction is that of 
the time of appearance for the peak of the response. 

The result of the simulated response for N=7 Is 
given in Figure 23* The plot shows several large amplitude 
oscillations and a peak height which is approximately 20 per- 
cent lower than the exact response. However* as compared to 
the N=3 case, the improvement is substantial. As far as 
computer time, the use of approximately 25 minutes on the 
IBM 360/50 was not too acceptable. This run formed a basis 
for a later comparative run on the CDC ?600, As table 8 
shows, the gain in execution speed with the CDC machine for 
N=7 Is approximately 75 times* 

The result of the simulated response for N=15 is 
given in Figure 24. The plot shows an increasingly good 
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agroement Kith the exact response* The oscillations are still 
present but of much reduced amplitude* The peak height Is 
slightly smaller than the exact height* However, the behavior 
of the response for 6<25 and 0>6o should be improved. This 
result indicated that a run with an additional number of col- 
location points was desirable. However, the large amount of 
computer time expended for the N=15 simulation, 110 minutes, 
was a debit on the side of further simulation* A small time 
run (10 minutes of computer time) for the N~21 case bn the 
IBM 360/50 gave an extrapolated complete run time (integrating 
up to 90 units of dimensionless time) of approximately 390 
minutes. This small time run indicated that the higher order 
approximation was decreasing the amplitude and frequency of 
initial oscillations as compared to the N=15 rxm* However, 
the time required to perform the complete calculations was 
too long (cost and scheduling) to obtain results on the 
IBM 360/50. 

At this point in time, effort was directed to ob- 
taining access to a computer more suited to the type of com- 
putations being made. Arrangements were made to remotely 
access the CDC 76OO computer at Combustion Engineering In 
Windsor, Connecticut, This machine's capabilities yielded a 
radical Improvement in expended computer time. The case of 
N«21 was run using' this machine. As Table 8 shows the run 
time to be approximately 122 times faster than the estimated 
run time for the reduced output case. Figure 25 shows the 
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results for N=21 . The results, when superimposed on the exact 

I 

response, give virtual exact agreement* The ojnly Identifiable 
discrepancies are the oscillations for 6<20 arid 0>?O* The 

i 

largest magnitude of the noted oscillations isj quite small, 

0.08 X icr^. 


Thus, it appears that for accurate approximation 
solutions for this problem, a fairly high degree of orthogonal 
collocation discretization in the spatial domain is required. 
Also, it is apparent that available computer hardware must be 
a very sizable consideration for extensive simulations* 

As a postscript to this part, it should be noted that 


going from the IBM 360 / 5 O to the CBC 76 OO required some alter- 
ations in the IMSL subroutine DBEBS, as the IMSL - CDC package 
was) not available at the Combustion Engineering CDC ?600, 

One of the changes Involved a machine-precision dependent con- 
stant. Fortunately, a CDC listing of DREBS was obtained in- 
directly from IMSL (Larsen, 19V^)* 



PART 8 

APPLICATION OF ORTHOGONAL COLLOCATION TO AN 
INTER-INTRAPARTICLE ADSORPTION MODEL V/ITlij 
NEGLIGIBLE MASS TRANSFER RESISTANCE j 

\ I 

i 

Orthogonal collocation approximations for the solu- 
tions of the previous two models have used discretizations in 
the interparticle region or axial dimension. When the model 
Is on© where concentration gradients are assumed to exist 
within the intrapartlcl© region, an approximation treatment 
for the Intraparticle domain is necessary. If one considers 
tho Inter-intrapartlole model represented by equations (2), 
(3)» (^)t a^id (5)» and considers tho case where N^q(j , 
the follov/lng model may ba deduced (see Appendix A)j 

(l/P0g)^>V/^z^ ~ ^y/dz - [(3(I-^)PA)(X4/r)Vpo^] = 

dy/dQ (92) 

(L/H)V?e 

...(93) 

(l/Rj)J^x^/6e = (9^) 

y* = m (95) 

This model is one of the most complicated forms 
that one might encounter In isothermal, packed bod analysts. 
The initial and boundary conditions are the same as equations 
(6) through (12) with equation (13) replacing (10) as the 
appropriate inter-lntra particle boundary condition. This 
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modification In the original model gives the inter- Intrapar- 
tide adsorption model negligible mass transfer resistance be- 
tween the inter particle and Intraparticle regions* For pur- 
poses of numerical treatment, the terminal boundary condition, 
equation (12), is replaced by the finite boundary condition 
previously applied. The forcing function, <J(0), is replaced 
by a finite width pulse, 0(e). 

This model introduces the necessity to perform a 
collocation approximation in the radial (intraparticle) di- 
mension, yu addition to the axial (Interparticle) dimension, 
z. For purposes of such a treatment, the trial function used 
by others (Finlayson, 1972, p*99) in the analysis of unsteady 
diffusion in a sphere, is equation (?0): 

Na 

= h(e) + (l-^)^^a3^(e) 

where Na is number of Intraparticle collocation points. 

When used in combination with an axial treatment, the axial 

position, z, should be Included to give: 

Na 

y^(z,^,o) = h(z,e) + (1-/^?) £ aj^(z,e) (96) 

^ 1=1 


The polynomials in JJ?" in equation ( 96 ) can be defined by a 
condition similar to equation (45) by making the change in 
variable /T^=z and 2ad/t=dz. The result is; 






(97) 
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In general, w(z) In equation (^5) Is of the formt 
w(z) = z^(l-z)^ 

Substituting this and the above change of variable Into 

equation (97) gives: 

b 

j = (01/2)^^ (98) 

The form given by Vllladsen ( 1970 ), The forjBulas used from 
that text for recursive computation of the respective poly- 
nomial coefficients defined by eqviatlon (45) may be modified, 
as Vllladsen shovrs, to give formulas for recursive computation 
of the coefficients for the polynomials in the squared dimen- 
sion defined by equation (98)* In the analysis that follows, 
the polynomials used are those defined by equation (98) with 
as=0, b=l,^=J,^ andj^=l. This is the case for spherical sym- 
metry, The coefficients are computed using the recursive 
formulas due to Vllladsen and the roots are taken from the 
values reported by Flnlayson ( 1972 , p, 102 ), 

Solution of this three, coupled partial differential 
equation problem requires orthogonal collocation approximations 
in two different spatial domains. The problem is one with 
three independent variables, z, A, and 0; and three dependent 
variables y, y , and yj (or x )♦ A solution to this type of 
problem appears not to have been attempted using the orthogonal 
collocation technique. To aid in envisioning the two-domain 
discretization required in the analysis of this problem. 
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Figure 26 gives a representation of the orthogonal collocation 
”grld”« Hero the collocation points are denoted by % for the 
Interparticle domain and for the intraparticle domain, 
respectively* 

Appendix J presents the orthogonal collocation ap- 
proximation development for the Intraparticle region. A com- 
puter program listing which was used to generate the first and 
second derivative Intrapartlole approximation matrices also 
appears- In this appendix. Appendix K develops expressions for 
simulating transient diffusion and adsorption/desorption be- 
havior within a single particle. Appendix L combines these 
results with the appropriate interparticle results to give the 
full representation of the orthogonal collocation approxima- 
tion for the model considered in this part. The result is a 
system of (Ne + 2 (Ke)(Wa)) coupled, ordinary differential 
equations: 

(He X 1) y a W (y - lg(e)) - COUPLE (99) 

vector — : 

and for ^=1 t . . ♦ t 

■ (Ha X 1) j = INTRA , - IHTRAC (1 y( 2 j. 6 )) 

vector ’ * 

+ IHTilAA y* ■ . (100) 

(Ha X 1) i? ^ = INTOM (x, , - x^ ,) (101) 

vector 

The strategy for determining what degree of approx- 
imation is adequate for accurate model solution is different 
than what was previously used; i.e . , comparislon of approximate 
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solutions to the exact solution. Here, successive approxima- 
tions must be compared to see If a convergent trend is noted. 
The appropriate order If approximation is then determined 
when increasing order gives no noticeable change in simulated 
response. This was the strategy which was to be used for an- 
alysis of this problem. However, as will be shown, completion 
of this strategy was not feasible,. 

The data set that is to be used in this section is 
the same that was used In the previous part for the rate of 
adsorption limited model. The differences are in the two 
choices of parameters mR^ and The parameter mHj was 

chosen to be based on an ioEq of 0 . 029 , a bed void fraction, 

€, of 0.40 and a particle void fraction, yS , of 0,40. Using 
equation ( 1 ), this would give the value of raR^ to be 0,01?4. 
The value of was estimated from Figure 19.^0 be 145.0 
from the model curve with a finite How^ever, because 

of the high in this data set, its contribution to the 

model variance is quite negligible and hence this value, of 
%U appropriate for the case of 

Table 9 presents a summary of what combination of 
interpartlole and Intraparticle approximations were slated for 
simulation. With access to the CEC ?600, it was decided to. 
run small-time (10 minutes of computer time) simulations on 
the IBM 360/50 to gain an estimate of the computer time neces- 
sary to complete the planned analysis, A listing of the pro- 
gram used for this purpose is given in Appendix L. Again 
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TABLE 9 


Inter-Intrapartlcle Model Orthogonal Collocation 
Approximations - Computer Time Estimates j 


Ne 

^A 

N 

Estimated Execution 
Time (min* )* 

3 

1 

9 

630.0 

3 

3 

21 

22500.0 

7 

1 

21 

1442.0 

7 

3 

49 

90000.0 

15 

- 1 

45 

3750.0 


K = ^2(1 + 2N^) 

» IBM 360 / 50 , F0RTRAN H, Integration (single precision) 
up to 90 units of dimensionless time v^ith output at 
approximately every 0*01 time units. 
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the IMSL subroutine CREBS was used to perform the Integration 

I 

of the system of ordinary differential equations for this 
problem. The time estimates in Table 9 ure overwhelming even 

I 

If one decreases them by the gain in speed (on; the order of 
one hundred times) with the CDC ’^ 600 ^ Within the imposed flnan- 
Dial constraints and justifications needed to motivate such 
V an expenditure, the complete analysis of this problem was 
not feasible. One could have possibly improved the situation 
by choosing “nice” parameter values but this would have 
negated the objective to study a technique with real-life 
problem conditions* 

The question arises then as to what information can 
be gained from this part of this investigation. First of all, 
it must be said that based on the computer hardware available 
and the technique used, the straightforward analysis of this 
complex problem using orthogonal collocation is not very ex- 
pedient. The pulsed, distributed system with multi-coupled 
transport phenomena presents a complicated problem for analysis. 
However, one might conceivably apply successfully the two do- 
main collocation approximation treatment for steady state or 
step response simulations for packed bed systems, either iso- 
thermal or non- isothermal. In addition, this analysis was 
based on two specific polynomial sets each orthogonal over 
one spatial domain interval. It could be argued that perhaps 
polynomials orthogonal to two domains simultaneously; * 

a surface, would be more appropriate for this type of problem. 
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Although the results of this section are in the 
negative side, they should not he construed as a condemnation 
of orthogonal collocation as a technique but rather as an 
example where a method may not be as well suited as others 
and where some innovations in the theory could possibly re- 
sult In a treatment that may be practical within the computa- 
s tlonal capabilities that now exist* 



PART 9 
DISCUSSION 

The initial part of thle investigation was motivated 
by previous efforts in the area of mathematical modeling of a 
gas chromatogi-aph. This Investigation set out to derive and 
study a model which Incorporated more of the dynamically rel- 
evant transport phenomena thought to be occurring in the ex- 
perimental systems being studied by Mars project co-v/orkers. 
Hence, a model has been proposed and derived which includes 
both interparticle and intraparticle transport phenomena. 

As with previous models, it was desirable to deter- 
mine whether or not this model could be used to adequately 
predict chromatograph system responses. Application of Laplace 
transform techniques gave a transform which was not readily 
invertable. However, because the model was linear and trans- 
formable, the derived transform could be used to determine 
the predictive capabilities of the model in the time domain. 
Here, the techniques of mo3nent analysis were applied arid it 
was shown that the model possessed a high degree of flexibility 
in predictive capabilities using the statistical quantities 
knovrn as moments which can characterize the responses of 
pulsed, distributed systems. This method of analysis is very 
useful because the effect of varying system parameters present 
in the model can be studied very efficiently and a great deal 
of insight into the model characteristics can be gained, as 
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was the case in this investigation* In fact, the results of 
the moment analysis gave sufficient motivation for the contin- 
ued analysis of the complex inter-lhtraparticle adsorption 
model. 

Because the derived model appeared to have no 
direct analytical solution, it V7as necessary to develop nu- 
merical capabilities in order to efficiently simulate the 
time domain response of the model and hence verify the model’s 
predictive effectiveness. However, prior to the investigation 
of numerical techniques, some study was given to the problem 
of replacing the infinite column boundary condition used in 
analytical work with a finite column terminal boundary condi- 
tion used in analysis of chemical reactor problems and which 
was necessary for numerical treatment of the model partial 
differential equations. It was desirable to apply the finite 
column boundary condition so that infinite column behavior 
would result at the bed outlet. Two simple problems were 
studied, one with the infinite column condition and the other 
with the finite oolvunn condition. Again, the technique of 
moment analysis proved a very effective tool in determlng how 
Infinite column response characteristics (moments) at the bed 
outlet might be matched by the problem with a finite terminal 
boundary condition. For the simple model considered, criteria 
were developed as a function of the Peclet number which gave 
guidelines for applying the finite terminal boundary. These 
gave reasonable assurance that column responses for both the 



finite and Infinite boundary conditions were essentially the 


same. These criteria were used, somewhat con|servatlvely, in 

the numerical solutions which were later made,! Although the 

i 

method of developing the boundary condition criteria was ap- 
plied to the simple model with one system parameter, the Peclet 
number, the method of analysis could be extended to more com- 
plex linear models and sets of criteria could be developed as 
functions of the system parameters present In each Individual 
model. 


Based on the preceding aspects of this .Investiga- 
tion,' effort was directed to the study of orthogonal collo- 
cation as a numerical approximation technique which would 
hopefully prove useful as an efficient tool for routine anal- 
ysis of the complex chromatograph system models. These 
models might be linear (as was the case in this work) or non- 
linear partial differential equations. The study of non-linear 
composition effects is an area of interest for continued 
chromatogrrph modeling effort. In this investigation, orthog- 
onal collocation was applied to approximate solutions to 
three linear, distributed model of Increasing complexity. 

The first model was a simple, one equation model requiring a 
collocation treatment in one spatial domain, the Interparticle, 
The second model Involved solution of two coupled partial 
differential equations requiring a collocation treatment In 
the Interparticle domain. The third model involved solution 
of three coupled partial differential equations requiring 
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collocation treatments In two spatial domains t the Interpar- 
tlcle and Intraparticle* 

The procedure for attacking these mathematical 
problems is summarized in Figure 27* This schematic provides 
a unified framework for discussing the general aspects and 
specific results of this investigation Into the use of orthog- 
onal collocation* 

The block denoted by STEP 1 serves as a starting 
point in problem analysis* This block, as indicated by the 
dashed lines is preliminary in nature and central to any mod- 
eling problem* For this Investigation, the work of Part 2 
could be lumped into this category. 

The block denoted by STEP 2 corresponds to that 
part of an analysis where one has to make a choice of the nu- 
merical method (If required) to use in the analysis of the 
formulated problem(s). The choices could be a finite differ- 
ence treatment, a finite element treatment, a weighted resid- 
ual treatment orthogonal collocation), or a varlation^il 

treatment. This choice may be motivated by previous exper- 

m 

ienoe, the work of other investigators in solving similar pro- 
blems, and/or the desire to establish the applicability of a 
certain method to a certain type of problem. In this inves- 
tigation, the choice of orthogonal collocation as a method of 
analysis was motivated by all of the above - the inefficiency 
of the finite difference technique to the simple, diffusion- 
convection problem (previous experience), the use of orthogonal 
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collocation by other investigators to solve certain problems 
in chemical reaction engineering, and the desire to ascertain 
the merits of orthogonal collocation as a computational tool 
for analysis of pulsed, distributed systems; , the gas 

chromatograph, ■ 

Following the choice of orthogonal collocation as 
the method for the analysis of the formulated problems, one 
proceeds to STEP 3, the choice of the trial function. Inputs 
to this choice may be innovation or imagination, suggestions 
from similar problems vjlth analytic solutions, or trial func- 
tions from previously worked ezaicples. The trial functions 
used in this investigation were taken from the work of other 
investigators. Hoviever, the trial function for interparticle 
analysis was generalized to an extent (Appendix G). The 
generality of this trial function was not explored - thorough 
exploration of the effects of trial function choice in com- 
bination with orthogonal function choice (STEP 4) for even 
one problem would be a basis for an entire investigation at 
least. 

This brings one to STEP 4, the choice of the or- 
thogonal fimctlons to be used in the trial function expansion. 
This block in the problem analysis can have the highest degree 
of variation* The choice can be dictated by the trial func- 
tion itself, symmetry considerations (the polynomials in 
for the intraparticle region), the type of solution (perhaps 
suggested by physical reasoning ), previous experience 
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(orthogonal polynomials uelghted In a certain way gave faster 
convergence with a previous problem), and the work of other 
Investigators, This investigator chose very specific poly- 
nomials for use In the trial functions employed in the pro- 
blems analysed. The choice was suggested by the iforks of 
other investigators and was further dictated by the availabil- 
ity of the required polynomial roots. Use of some less 
** common** polynomials require determination of roots which may 
not be tabulated to high accuracy. Thus, one would have to 
pursue root-finding computations vdiloh, based on available 
computer precision capabilities, may be infeasible. This 
type of study; i*e . , different polynomial types, was a desired 
component of this investigation but efforts to compute nevf 
roots to sixteen figure accuracy were limited by the available 
IBM 360/50, Thus, this desired area of study vras abandoned. 

In conjunction with this, one should note the added input to 
the ST:ip 4 block citing the very practical consideration of 
computing capabilities - in this Instance, woi'd-length capa- 
bility, Recently, the notion of there being better polynom- 
ials for certain problems received attention by Ramkrlshna 
(1973)0 He showed that the choice of “problem specific poly- 
nomials’* to be relevant and desirable for effective use of 
weighted residual techniques. 

The block denoted by STEP 5 is labeled DISCRETIZATION. 
This is descriptive of the manipulations and computations 
which must be made to reduce the expressions for the partial 



109 


derivatives at the collocation points to expressions in terms 
of the solutions at the respective collocation points. The 
manipulations of Appendix G and Appendix J are typical of 
what must be done, ' A key aspect of these computations Is the 
inversion of a matrix. Part 5 showed how the conditioning of 
the matrix to be inverted may be enhanced by a change In for- 
mulation. One was able to retain a tolerable condition with- 
in the constraint of the available computer precision. How- 
ever, as v/as pointed out in Part 6, attempts to generate a 
thirty-first order interparticle discretization were unsuc- 
cessful because the word-length of the available IBM 3^0/50 
computer limited the accuracy of the polynomial coefficients. 

Next is STEP SIMULATION. This block comprises 
the use of the previously derived and computed discrete re- 
presentations to reduce the distributed model to a set of 
ordinary differential equations. This set of ordinary dif- 
ferential equations can be integrated to yield the approx- 
imate response. As was done with the simple model (Part 6), 
the equations can be put In a suitable form. wherein an eigen- 
analysis of the system can be made to determine the character 
of the approximation solution* This also served to expose 
the stiffness of the equation set. The simulated response(s) 
can be compared v^lth exact solutions (if available), solutions 
from other techniques tif available), and with simulations 
using different orders of approximation. As was shown with 
the rate of adsorption limited Inter-intraparticle model, the 
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available computing capabilities (execution speed) can be a 
factor In limiting the extent of any planned simulation pro- 
gram. In fact, this constraint (even with the CCC ?600) 
prohibited complete analysis of the Inter-intrapartlcle ad- 
sorption model with negligible mass transfer resistance be- 
tween the interparticle and Intraparticle regions. 

^ STEP 7, ASSESSMENT, serves as an area where one 

can assess the results and the reasons for the results. Among 
the questions that must be answered are; 

1 . Is the problem practically solvable? 

2. Does the orthogonal collocation treatment, as 
applied , have sufficient power to be used as a 
routine tool In analysis of the posed problem(s)? 

Regarding the first two problems solved in this’ investigation, 
the answers to the above questions would be affirmative. How- 
ever, with regard to the third problem, the answers are not 
affirmative. The key words In the second question are “ as 
applied ** because the trial function and/or polynomial type 
may be unsulted to the problem at hand and may thus require 
some new Innovations In this area. This is the reason for 
the “feedback** loops from STEP 7 to STEP 3 stnd to STEP 4. 
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CONCLUSIONS 

This Investigation has been conducted in conjunction 
Kith the group effort to define fundamental system design 
criteria necessary for an optimal design of a combination gas 
chromatograph - mass spectrometer* Specifically, this inves~ 
tlgatlon has dealt with the formulation of a more complex 
mathematical model for a gas chromatograph and subsequent ef- 
forts to ascertain the merits of the numerical technique knovjn 
as orthogonal collocation as a technique v^orthy of I'outine use 
In the time domain simulation of complex gas chromatograph 
models* 

Previous work dictated the formulation of a model 
which took into account more of the dynamically relevant 
transport mechanisms thought to be occurring In the chromato- 
graph system^ A model has been formulated v/hich Includes 
Intraparticle diffusion and rates of adsorption that were 
heretofore neglected. The model has been analysed using the 
moment analysis technique. This analysis of the proposed 
Inter-Intraparticle Adsorption Model indicates that the gross 
characteristics of actual data are more adequately predicted 
than with previous models. 

The mathematical complexity of the proposed Inter- 
Intrapartlcle Adsorption Model has prompted consideration of 
numerical techniques appropriate for the solution of the 
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partial differential equations being postulated. The use 
of numerical techniques for the second-order models being 
considered requires the use of a finite terminal boundary 
condition. Criteria have been developed for a simple model 
wherein a finite terminal boundary condition can be applied 
which yields system responses which are for all intents and 
purposes equivalent to the responses obtained using an in- 
finite column boundary condition, ■ 

The general theory and computational aspects of 
the method loiovjn as orthogonal collocation have been reviewed 
and discussed. An alternate method of problem formulation 
gives a matrix (vihich must be inverted in either formulation) 
which is significantly better conditioned for Inversion pur- 
poses. It is concluded that this different approach is better 
than previously documented approaches when computer word- 
length capabilities are a consideration as Is the case for 
most practical situations. 

The method of orthogonal -collocation has been suc- 
cessfully applied to two problems of the chromatograph system 
type. The first problem was the simple transient diffusion- 
convection equation and the second was the rate of adsorption 
limited Inter-lntrapartlcle model. These models required 
orthogonal collocation treatments for one spatial domain, the 
Interpartlcle, For the system parameters considered it appears 
that 15 collocation points are adequate for the simple model 
and 21 collocation points are adequate for the rate of ad- 
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sorption limited model. However, the latter model required 
the computing power of a CDC 76 OO. 

The application of orthogonal collocation to an 
Inter-lntrapartlcle adsorption model with negligible mass 

t 

transfer resistance between the interparticle and Intrapar- 
tide regions is not, based on the two domain (interparticle 
and Intraparticle) orthogonal collocation treatment, prac- 
tically feasible even with the computing power of a CDC 7^00* 

It Is concluded that although the theory of orthogonal col- 
location may be viable, there could possibly be significant 
improvement in practical requirements if modifications in 
trial function and/or orthogonal function choices can suc- 
cessfully be effected. This conclusion applies, to a lesser 
degree, to the two other models considered in this investi- 
gation. 

Thus, It appears that in its present state, orthog- 
onal collocation can be a useful tool for analysis of one 
spatial domain, pulsed, distributed systems. Use of orthogc lal 
collocation for two-spatial domain, pulsed, distributed 
systems requires the reversion back to the steps of trial 
function selection and orthogonal function selection in order 
to effect a practical approximation treatment* 

Throughout this investigation certain areas of work 
have been mentioned, as areas suitable for future research. 

The proceeding remarks summarize these areas. 

One area is the use of specific polynomials to take 





advantage of their orthogonality relative to the trial func- 
tion expansion and whether specific polynomials within the 
trial function would produce better results from the stand- 
point of Increased accuracy with a lower order of approximation* 
In addition,, the form of the trial function is an area where 
further investigation may be made to ascertain what trial 
.function form (in conjunction with orthogonal function choice) 
is best for a given problem. 

The solution of the sets of ordinary differential 
equations produced by application of the orthogonal collo- 
cation technique is another area suitable for further re- 
search, It vfas shown that the equations for the orthogonal 
collocation approximation of the simple, diffusion-convection 
model possess characteristics of a stiff set* This situ- 
ation raises the question as to what method of integration 
should be used. This could form an additional area of re- 
search - whether orthogonal collocation approximations pro- 
duce, in general, stiff ordinary differential equation sets 
and what integration algorithm can be used most effectively 
for simulation purposes. 
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NOMENCLATURE 

- unit impulse j Dirac delta function. 

- areas under output response curve, input 
response curve, and impulse curve, respec- 
tively. 

- lower bound on interval of orthogonality 
used in orthogonal polynomial definition, 
equation (45) • 

- time-varying coefficients in trial function 
expansion. 

« ratio bf interfaolal area to packed volume. 

- upper bound of interval of orthogonality 
used in orthogonal polynomial definition,, 
equation (45)* 

- scale factor used in orthogonal polynomial 
definition, equation (45). 

vector used in equation (99). 

- time-varying coefficients in trial function 
expansion* 

vector of time-varying coefficients defined 
In equation ( 63 ). 

- time-varying function in the trial function 
expansion. 

- vector of time-varying coefficients defined 
in equation ( 50 ). 
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g(e) 

h(0) 


IHTRA 

IMTRAA 

INTRAC 

IMTRAE 
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m 
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% 

NhU 

^tOG 

Pe 
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- time- varying function in the, trial function 

I 

expansion, j 

- time-varying function in equation (?0), Later 

i 

extended to h(z,0) in equation ( 96 ), 

- modified Bessel function of the first kind, 

- matrices used in equation (100), 


- matrix used in equation (101), 

- length of chromatograph column^ dimensional, 

- r equilibrium constant, 

- number of collocation points except as defined 
differently in Part 8. 

- number of Interparticls collocation points, 

- number of intraparticle collocation points, 

- the number of reactor units, a dimensionless 
measure of the rate of adsorption, 

- number of transfer units, dimensionless, 

- Peclet number, dimensionless, 

- intraparticle Peclet number, a dimensionless 
measure of diffusion rates vrlthin the particle, 

- Interparticle Peclet number^ a dimensionless 
measure of diffusion rates within the carrier 
gas. 

- group of polynomials, initially arbitrary 
but later constrained to be orthogonal on 
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interval ra,b1 by equation i(45) or (98)* 

■ , . ! 

- matrix defined by equation (| 62 ). 

- matrix defined by equation (!66). 

i 

- matrix defined by equation ('67). 

- Intraparticle space variable, dimensionless* 

- particle radius, dimensional* 

- rate of sample adsorption. 

- moles of fluid In the particle per mole of 
adsorption sites v;lthin the particle* 

- moles of fluid within the total bed per moles 
of adsorption sites within the bed* 

- matrix defined by equation (49)* 

- matrix defined by equation (5^)« 

- matrix defined by equation (57)* 

- residual formed by trial function substi- 
tution in a differential equation* 

- I^place transform variable* 

- weighting function used in orthogonal, poly- 
nomial defining equation (45). 

- weighting function in weighted residual 
Integral, equation (43). 

- matrix used initially in equation (82); then 
in equation (90). 

- matrices used in equation ( 90 ) 


- matrix used in equation (91). 
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y 


yi 


n 


1 % 




z 

2 ;) 


adsorbed phase concentration, dimensionless. 
Interparticle gas phase composition, dimen- 
sionless, 

Intraparticle gas phase composition, dimen- 
sionless, i 

equilibrium intraparticle gas phase compo- 
sition, dimensionless* 

vector of solution values at the Interparticle 
collocation points* 

vector of equilibrium concentration values 
at the Interparticle collocation points* 
vector of intraparticle concentration values 
at the intraparticle collocation points* 
vector of intraparticle equilibrium con- 
centration values at the intraparticle 
collocation points. . , 

axial position in column, dimensionless* 
collocation point or end point, dimension* 
less* 

axial position vjhere finite terminal boundary 
condition, equation (33)* is applied* 


GREEK LETTERS 




» .. part of the power of 




1 


in equation (98); 



particle porosity or void fraction; power of 
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AQ 

Az 

6 ( ) 
6iJ 

7T 


Mn 


0(e) 


quantity (1-^^) In equation (98) • 

- tlmo increment In finite difference method* 

- space Increment in finite difference method. 

- Dirac delta function. 

- Kroneoker delta. 

- vttld fraction of the "bed. 

- 3 * 141 59 ... 

- dimensionless time variable. 

th 

« the h moment about the origin defined by 
equation (14). 

- the n moment about defined by 

equation (1?). 

- function which satisfied boundary condition 
In trial function expansion. 

- forcing function used In analysis of chro- 
matograph problems. 


■ SUBSCRIPTS 

I - refers to Case I boundary condition analysis. 

II - refers to Case II boundary condition analysis. 

1 - refers to space level in Finite Difference 

technique; refers to column in Orthogonal 
Collocation matrices. 

j - refers to time level in Finite Difference 

techniques; refers to row and/or collocation 
points in. Orthogonal Collocation matrices. 
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MISCELLA.WEOUS 



« refers to the matrix element of the 

+*Vi 

and the column* 


row 
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